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A.  SCIENTIFIC  OBJECTIVES:  The  objective  of  this  research  program  is  to  apply 
model-based  signal  processing  methods  to  enhance  the  performance  of  computation 
electromagnetics  (CEM)  simulators  for  shipboard  antenna  design.  While  recent  advances 
in  CEM  algorithms  has  significantly  reduced  the  simulation  cost  of  modeling  complex 
radiation  and  scattering  phenomena,  real-world  engineering  design  and  optimization 
often  require  that  calculations  be  carried  out  repeatedly  over  large  parameter  spaces  such 
as  frequency  and  aspect  angle,  making  the  computation  cost  still  exceedingly  high.  The 
goal  of  this  research  is  to  apply  model-based  signal  processing  algorithms  to  CEM 
simulators  to  overcome  such  computational  bottleneck  and  achieve  higher  design 
throughput.  In  particular,  we  address  the  issue  of  how  to  interpolate  and  extrapolate  the 
frequency  and  angular  behaviors  of  antenna  characteristics  based  on  a  sparse  set  of 
computed  data.  We  develop  algorithms  to  extract  sparse  and  physical  models  of  the 
relevant  physics  embedded  in  CEM  data  to  facilitate  design  and  synthesis.  Finally,  we 
explore  the  transportability  of  the  methodology  to  other  applications  such  as  radar 
signatures  and  wireless  channel  characteristics. 

B.  SUMMARY  OF  RESULTS  AND  SIGNIFICANT  ACCOMPLISHMENTS:  Our 
research  during  the  three-year  duration  of  this  program  has  been  focused  on  three  topics. 
First,  we  have  developed  a  number  of  model-based  algorithms  to  extrapolate  and 
interpolate  the  frequency  and  angular  behaviors  of  antenna  radiation  and  scattering 
characteristics  on  complex  platforms.  These  algorithms  can  lead  to  orders  of  magnitude 


reduction  in  computation  time  in  populating  a  large  data  set  from  a  limited  number  of 
computed  data.  Second,  we  have  developed  algorithms  to  extract  sparse  and  highly 
physical  models  of  antenna  radiation  data  on  complex  platforms  from  CEM  data.  The 
extracted  features  can  be  used  to  pinpoint  the  locations  of  platform  scattering  for 
engineering  design  and  optimization.  Third,  we  have  investigated  the  effect  of  mutual 
coupling  and  platform  interaction  on  the  performance  of  antenna  arrays  mounted  on 
complex  platforms.  Both  theoretical  analysis  and  measurements  have  been  earned  out  to 
assess  the  near-field  effect  of  the  mounting  environment  on  the  direction  finding 
performance  of  antenna  arrays.  In  addition  to  these  three  main  topics,  we  have  also 
carried  out  exploratory  research  in  two  areas.  First,  we  have  proposed  and  implemented 
algorithms  to  improve  the  speed  and  accuracy  of  CEM  solvers  based  on  the  fast  multipole 
method.  We  have  proposed  a  new  set  of  triangular  basis  functions  to  improve  the 
accuracy  of  the  traditional  Rao-Wilton-Glisson  basis.  We  have  also  proposed  a  low- 
complexity,  wavelet-based  preconditioner  to  accelerate  the  convergence  rate  of  iterative 
solvers.  Second,  we  have  explored  the  use  of  genetic  algorithm  together  with  CEM 
solvers  to  design  antenna  structures  for  optimal  performance.  This  work  leverages  the 
results  of  this  research  program  for  actual  design  synthesis.  Below  we  shall  describe  the 
five  topics  in  more  detail.  Emphasis  will  be  placed  on  our  progress  during  the  last  year  of 
the  program. 

Model-Based  Extrapolation  and  Interpolation  of  Radiation  and  Scattering  Data 
for  Complex  Platforms.  In  antenna  design  and  analysis,  the  mounting  platform  can 
have  significant  impact  on  the  resulting  antenna  radiation  characteristics.  Rigorous 
numerical  solution  of  the  radiation  problem  over  a  complex  platform  is  very  time 
consuming,  and  the  computation  complexity  increases  dramatically  as  the  frequency 
increases.  During  the  first  year  of  this  program,  we  developed  a  model-based  frequency 
extrapolation  technique,  with  which  the  radiated  field  over  a  broad  band  of  frequencies 
can  be  obtained  using  rigorously  computed  results  at  low  frequencies.  Our  approach 
entailed  fitting  the  currents  computed  at  low  frequencies  to  a  time-of-amval  model  and 
determining  the  model  coefficients  using  the  superresolution  algorithm  ESPRIT  [1].  The 
currents  and  the  radiation  characteristics  were  then  extrapolated  from  the  resulting  model. 


Some  initial  results  were  obtained  during  the  first  year  but  the  accuracy  of  the 
extrapolated  results  was  not  satisfactory. 

During  the  second  year,  we  improved  the  algorithm  significantly  by  adopting  a 
frequency-dependent  model  for  the  induced  current.  The  different  time-of-amval  terms 
correspond  to  different  incident  wave  mechanisms  from  both  the  direct  antenna  radiation 
and  higher-order  scattering  from  other  parts  of  the  platform.  These  different  mechanisms 
in  general  have  different  frequency  dependencies.  For  canonical  platform  shapes,  their 
exact  frequency  dependencies  are  known  through  the  geometrical  theory  of  diffraction 
(GTD)  [2].  However,  for  more  complex  structures,  they  must  be  determined  numerically. 
We  devised  a  pre-multiplication  scheme  in  conjunction  with  the  complex  time-of-amval 
estimation  from  ESPRIT  to  determine  the  additional  frequency-dependent  factors  [27]. 
The  performance  of  the  algorithm  in  the  presence  of  noise  was  evaluated  based  on 
simulated  data  and  errors  in  the  estimation  of  model  parameters  were  quantified.  Our 
results  showed  that  the  method  was  quite  robust.  The  algorithm  was  also  extended  in 
order  to  simultaneously  extrapolate  frequency  and  aspect  data  in  scattering  problems 
[19]. 

We  also  investigated  the  interpolation  problem  for  CEM  data.  Although  the  potential 
payoff  of  an  interpolation  algorithm  is  not  as  great  as  extrapolation,  it  provides  a  more 
robust  way  to  achieve  computational  savings.  To  achieve  the  parameterization  in  the 
frequency  dimension,  we  again  utilized  the  multiple-arrival  model  for  the  induced 
current.  To  obtain  the  model  parameters  from  the  computed  data,  we  adopted  an 
algorithm  termed  adaptive  feature  extraction  (AFE).  This  algorithm  was  originally 
devised  by  us  to  construct  the  inverse  synthetic  aperture  radar  (ISAR)  image  from  radar 
measurement  data  that  was  undersampled  in  the  aspect  dimension  [15].  The  essential  idea 
of  the  AFE  algorithm  is  to  search  and  extract  out  individual  scattering  features  from  the 
data  set  one  at  a  time.  During  each  iteration,  the  strongest  feature  is  identified  and 
removed  from  the  original  data.  The  procedure  is  then  iterated  until  the  data  is  well 
parameterized  by  the  feature  set.  In  this  manner,  the  interference  between  different 
scattering  features,  which  is  significant  for  undersampled  data,  can  be  largely  avoided. 
Since  the  features  in  our  model  are  exponential  functions  of  frequency  and  aspect, 
random  sampling  was  used  during  the  collection  of  the  original  computation  data  in  order 


to  avoid  the  ambiguity  in  selecting  the  strongest  feature.  This  algorithm  was  extended  to 
interpolate  2-dimensional  frequency  and  angle  data  [26].  In  this  case,  a  more  general 
time-of-arrival  and  angle-of-arrival  model  was  used  to  parameterize  the  current.  The  2D 
algorithm  was  used  in  conjunction  with  the  fast  multipole  code  FISC  [3]  to  predict  the 
radar  image  of  the  benchmark  VFY-218  airplane  at  UHF  band  [4].  Excellent  agreement 
with  measured  data  was  obtained.  The  resulting  computation  time  savings  was  about  two 
orders  of  magnitude  lower  than  the  brute-force  computation. 

In  the  third  year,  we  focused  our  attention  on  the  resonant  behaviors  of  antenna- 
platform  configurations.  This  topic  was  motivated  by  the  interest  of  the  Navy  SPAWAR 
Center  in  HF  antennas  where  ship  body  resonances  could  dominate  antenna  frequency 
characteristics.  In  the  lower  frequency  regime,  the  scattering  phenomenology  differs  from 
ray-optical  characteristics.  Consequently,  the  time-of-arrival  model  we  had  utilized 
successfully  was  no  longer  efficient  in  this  regime.  We  proposed  a  hybridization  of  the 
multiple-arrival  model  and  the  rational  function  model  [21].  The  AFE  algorithm  was 
again  used  to  extract  the  parameters  in  the  multiple-arrival  model,  while  the  parameters  in 
the  rational  function  model  were  obtained  using  the  Cauchy  method  [5].  At  each  point  on 
the  target,  we  first  parameterized  the  induced  current  using  both  the  AFE  and  Cauchy 
algorithms  separately.  The  model  that  best  matched  the  scattering  physics  at  that  location 
was  automatically  chosen  based  on  the  resulting  interpolation  error.  Consequently,  those 
regions  on  the  target  that  could  best  be  described  by  ray-optical  phenomena  were 
interpolated  by  the  multiple-arrival  model,  while  those  regions  that  exhibited  resonance 
phenomena  were  interpolated  by  the  rational  function  model.  Numerical  results  for  a 
target  containing  complex  features  were  generated.  It  was  shown  that  the  hybrid  model 
resulted  in  more  accurate  interpolation  than  either  of  the  two  models  alone.  Extension  of 
this  algorithm  to  3-D,  non-conducting  targets  was  also  attempted. 

Extraction  of  Sparse  Representation  of  Antenna  Radiation  Data.  The  previous 
topic  is  focused  on  the  utilization  of  model-based  processing  to  extrapolate  or  interpolate 
CEM  data  for  the  “forward”  solution  of  electromagnetic  radiation  and  scattering 
problems  involving  complex  platforms.  An  equally  important  problem  from  the  design 
perspective  is  the  “inverse”  algorithm  of  spatially  pinpointing  the  locations  of  platform 


scattering  based  on  CEM  data.  During  the  first  year,  we  developed  a  method  to  extract  a 
sparse  model  of  the  antenna  radiation  pattern  on  a  complex  platform  [17].  This 
representation  is  based  on  a  point  radiator  model  that  describes  the  radiation  pattern  by  a 
collection  of  “radiation  centers”  on  the  platform.  The  methodology  for  obtaining  the 
radiation  center  model  entailed  first  generating  the  3-D  antenna  synthetic  aperture  radar 
(ASAR)  imagery  of  the  platform,  and  then  parameterizing  the  resulting  image  by  a 
collection  of  point  radiators  via  the  CLEAN  algorithm  [6].  Once  such  a  representation 
was  obtained,  we  could  rapidly  reconstruct  antenna  radiation  patterns  over  frequencies 
and  aspects  with  good  fidelity.  Thus  such  a  model  could  be  used  for  real-time 
reconstruction  of  complex  antenna  patterns  in  high-level  system  simulations. 
Furthermore,  the  resulting  radiation  center  information  could  be  used  to  pinpoint  cause- 
and-effect  in  platform  scattering  and  provide  important  design  guidelines  for  antenna 
placement  and  optimization.  Experimental  verification  of  the  ASAR  concept  was 
reported  recently  in  [7]. 

The  ASAR  method  is  a  Fourier-based  algorithm  that  relies  on  a  number  of  small- 
angle,  small-bandwidth  approximations.  Furthermore,  the  concept  was  only  demonstrated 
using  high-frequency  ray-tracing  simulation  but  not  more  rigorous  CEM  data.  During  the 
second  year,  we  overcame  the  above  deficiencies  by  developing  a  more  generalized 
matching  pursuit  algorithm  [8]  in  the  frequency-aspect  domain  based  on  a  point  radiator 
basis  [18].  We  also  demonstrated  the  feasibility  of  extracting  a  sparse  model  of  the 
antenna-platform  interaction  using  CEM  data  from  FISC.  The  matching  pursuit 
algorithm  was  implemented  based  on  the  radiation  center  basis.  To  speed  up  the 
parameterization  time  of  the  matching  pursuit  algorithm,  we  estimated  the  location  of  the 
radiation  centers  by  utilizing  the  Fourier-based  ASAR  algorithm.  The  point  with  the 
highest  intensity  was  first  located  in  the  ASAR  image  and  its  amplitude  and  coordinates 
were  determined  to  serve  as  an  estimate  of  the  strongest  radiation  center.  We  next 
zoomed  in  on  the  precise  location  of  the  radiation  center  via  a  fine  search.  We  then 
subtracted  the  contribution  of  this  radiation  center  from  the  total  radiated  field  and 
iterated  the  process  until  the  energy  in  the  residual  signal  had  reached  a  sufficiently  small 
level.  Our  results  showed  that  radiation  patterns  from  a  complex  ship-like  platform  could 
be  well  represented  by  a  very  sparse  set  of  radiation  centers. 


Finally,  we  also  investigated  ways  to  predict  the  near-field  radar  cross  section  (RCS) 
of  a  target  by  taking  advantage  of  the  far-field  scattering  centers  [14].  The  basic  idea  was 
to  extract  the  far-field  3D  scattering  centers  once  and  then  use  the  strengths  and  positions 
of  the  extracted  scattering  centers  to  reconstruct  the  near-field  RCS  at  different 
transceiver  locations.  Since  the  extraction  was  done  only  at  the  far-field  range  and  the 
reconstruction  could  be  carried  out  in  near  real-time  for  a  variety  of  sensor  parameters, 
this  method  is  a  very  efficient  way  to  generate  near-field  RCS  data.  Although  only  near¬ 
field  monostatic  scattering  was  considered,  the  method  could  be  easily  extended  to 
handle  more  general  bistatic  scenarios. 

Effect  of  Mutual  Coupling  and  Platform  Effects  on  Antenna  Arrays.  During  the 
third  year  of  this  program,  we  initiated  an  investigation  into  mutual  coupling  and 
platform  effects  for  antenna  arrays  mounted  on  complex  platforms.  First,  the  effect  of 
mutual  coupling  on  direction  finding  performance  of  circular  arrays  was  simulated  using 
electromagnetic  computation  and  the  MUSIC  superresolution  algorithm  [20].  The  various 
scenarios  studied  included  element  spacing,  noise  level,  antenna  gain  and  delay,  and 
platform  effects.  Our  results  showed  quantitatively  the  effects  of  mutual  coupling  on  the 
direction  of  arrival  estimates.  The  compensation  of  the  coupling  effect  using  the 
coupling  matrix  approach  was  also  examined.  This  approach  was  found  to  be  quite 
satisfactory  in  most  cases,  except  when  the  array  calibration  data  contained  high  noise 
level,  or  when  strong  platform  effects  were  present. 

As  a  follow-up  to  the  theoretical  study,  we  also  made  measurements  using  the  7- 
element  circular  array  testbed  operated  by  the  wireless  group  at  the  University  of  Texas 
[22].  Previous  experimental  studies  by  this  group  using  the  array  had  shown  anomalous 
DO  A  results  [9].  For  a  single  mobile  user  in  an  open  field  environment,  the  spatial 
spectra  using  the  array  showed  a  large  “rebounce”  sidelobe  not  corresponding  to  any 
signal  multipath  or  interference.  In  our  study,  we  computed  the  array  response  including 
mutual  coupling  effects  using  the  Numerical  Electromagnetics  Code  (NEC).  This 
computed  array  response  was  then  applied  to  DOA  estimation  using  measurement  data 
collected  from  the  antenna  array.  Our  results  showed  that  the  rebounce  sidelobe  was 
significantly  reduced  when  mutual  coupling  was  taken  into  consideration.  The  same 


study  of  two  co-channel  mobile  users  was  also  carried  out  and  improved  DOA 
performance  was  achieved  [29]. 

As  a  step  toward  investigating  platform  effects  on  array  performance,  we  examined 
two  commonly  used  approaches  to  determine  the  coupling  matrix  and  their  limitations 
[23].  The  coupling  matrix  obtained  from  the  minimum  mean  square  error  matching  of  the 
active  and  stand-alone  element  patterns  was  shown  to  be  more  effective  than  the  mutual 
impedance  approach.  Furthermore,  we  proposed  an  extended  coupling  matrix  model  for 
more  complex  antenna  arrays  where  parasitic  elements  were  present.  Finally,  we  devised 
a  combined  model  to  model  antenna  array  response  including  both  mutual  coupling  and 
platform  effects  [33].  The  mutual  coupling  was  described  by  a  coupling  matrix  while  the 
modification  to  the  incident  field  due  to  the  platform  was  modeled  using  the  point 
scatterer  model.  A  matching  pursuit  algorithm  was  employed  to  determine  the  model 
coefficients.  Numerical  results  showed  that  the  active  element  patterns  predicted  by  the 
model  were  in  good  agreement  with  the  rigorously  simulated  patterns.  The  advantage  of 
the  new  combined  model  is  that  it  provides  physical  insights  into  platform  effects  by 
pinpointing  the  locations  of  the  dominant  scattering  mechanisms  on  the  platform. 

Improving  the  Speed  and  Accuracy  of  Fast  Multipole  Solvers.  In  addition  to  the 
three  main  topics  described  above,  we  also  proposed  and  implemented  algorithms  to 
improve  the  speed  and  accuracy  of  CEM  solvers  based  on  the  fast  multipole  method. 
First,  we  devised  a  new  set  of  triangular  basis  functions  to  improve  the  accuracy  of  the 
traditional  Rao-Wilton-Glisson  (RWG)  basis  [10,30,32].  In  the  computation  of  surface 
currents  and  near  field  data,  we  found  that  the  tradition  RWG  triangular  patch  basis  did 
not  provide  sufficient  accuracy.  The  basic  underlying  problem  is  that  RWG  basis  is  not 
complete  and  thus  cannot  represent  an  arbitrary  linear  current  distribution  on  a  patch. 
Subsequently,  we  devised  an  enhanced  set  of  triangular  patch  basis  functions.  It  was 
shown  that  these  new  basis  functions  provided  a  much  better  representation  of  the  current 
distribution  than  the  RWG  basis  while  keeping  the  same  computational  complexity.  In 
the  examples  we  analyzed,  we  always  obtained,  for  the  same  triangulation  scheme,  a 
better  result  using  the  new  basis  functions  in  terms  of  RCS  prediction.  In  terms  of  current 
distribution,  the  comparative  performance  was  even  better  for  this  new  formulation.  For 


such  application  as  determining  the  near  field  distribution,  this  set  of  basis  functions  can 
result  in  significant  accuracy  improvement  over  the  traditional  RWG  basis.  The 
computational  complexity  involved  was  about  the  same  as  the  RWG  basis.  Therefore, 
existing  computer  codes  can  be  adapted  to  use  them  with  minimal  modifications. 

Second,  we  proposed  a  low-complexity,  wavelet-based  preconditioner  to  accelerate 
the  convergence  rate  of  iterative  CEM  solvers  [25,28,34].  As  pointed  out  earlier,  ship 
body  resonances  can  dominate  antenna  frequency  characteristics  in  the  lower  frequency 
regime.  An  important  problem  related  to  the  study  of  resonant  region  using  an  iterative 
CEM  solver  is  that  when  the  platform  exhibits  strong  resonance  effect,  the  iteration 
number  required  for  an  accurate  solution  can  become  large.  A  good  preconditioner  for 
the  moment  matrix  is  needed  to  alleviate  this  problem.  We  had  previously  investigated 
the  use  of  wavelet  packet  basis  for  the  sparsification  of  moment  matrix  [11,12].  We 
extended  our  work  by  devising  an  approximate-inverse  preconditioner  to  the  moment 
equation  using  the  wavelet  representation.  In  our  approach,  we  constructed  the 
approximate-inverse  preconditioner  in  the  wavelet  domain  where  both  the  moment  matrix 
and  its  inverse  exhibit  sparse,  multilevel  finger  structures.  The  inversion  was  carried  out 
as  a  Forbenius-norm  minimization  problem.  Numerical  results  on  a  3-D  cavity  showed 
that  the  iteration  numbers  were  significantly  reduced  with  the  preconditioned  system. 
More  importantly,  the  computational  cost  of  the  preconditioner  was  kept  under 
OiNlogN),  making  it  compatible  with  the  fast  multipole  method. 

Antenna  Design  and  Optimization  Using  Genetic  Algorithms.  As  an  exploratory 
effort  during  the  last  year  of  the  this  program,  we  initiated  a  study  into  applying  genetic 
algorithms  (GA)  for  the  design  and  optimization  of  antenna  structures.  Genetic 
algorithms  are  stochastic  search  methods  based  on  the  concept  of  natural  selection  and 
evolution.  This  class  of  algorithms  is  particularly  attractive  for  finding  an  approximate 
global  optimum  in  a  very  high-dimensional  space.  Research  in  the  application  of  GA  for 
antenna  design  has  been  ongoing  since  the  early  1990s  [13].  In  contrast  to  a  local 
optimization  algorithm,  GA  allows  the  design  space  to  be  more  fully  explored  (at  the 
expense  of  computational  cost).  In  the  GA  step,  the  parameters  to  be  optimized  are  first 
encoded  into  a  binary  string  known  as  a  chromosome.  Chromosomes  in  the  population 


are  ranked  according  to  a  chosen  cost  function  and  a  new  generation  of  children 
chromosomes  is  produced  through  the  rules  of  heredity  and  mutation.  This  process  is 
repeated  until  a  satisfactory  set  of  parameters  is  found  that  best  meets  the  design  criteria. 

We  carried  out  a  preliminary  GA  study  to  optimize  the  shape  of  microstrip  antenna 
for  broadband  operation  [24,31].  The  geometry  of  the  patch  was  discretized  into  a  two- 
dimensional  binary  map  and  GA  was  used  to  search  for  the  optimal  patch  shape  that 
maximized  bandwidth.  A  fast  CEM  code  specifically  tailored  for  microstrip  structures 
was  used.  Several  schemes  were  utilized  to  accelerate  the  convergence  of  the  GA 
including  2-point  crossover  and  geometrical  filtering  through  the  use  of  median  filter. 
The  optimized  shape  showed  a  four-fold  improvement  in  bandwidth  when  compared  to  a 
standard  square  microstrip  antenna.  This  result  was  verified  by  laboratory  measurement. 
The  basic  operating  principle  of  the  optimized  shape  was  interpreted  in  terms  of  a 
combination  of  dual-mode  operation  and  ragged  edge  shape. 

We  also  explored  the  use  of  GA  for  model  extraction.  In  particular,  we  addressed  the 
problem  of  finding  the  equivalent  impedance  boundary  condition  (IBC)  for  a  corrugated 
material  coating  structure  [16].  In  our  approach,  rigorous  solutions  of  the  reflection 
coefficients  at  a  number  of  incident  angles  were  first  calculated  using  a  periodic  method 
of  moments  (MoM)  solver.  The  IBC  model  was  used  to  predict  the  reflection  coefficients 
at  the  same  observation  angles.  The  model  coefficients  were  then  optimized  using  GA  so 
that  the  difference  between  the  approximated  and  the  MoM-predicted  reflection 
coefficients  was  minimized.  The  genetic  algorithm  proved  efficient  in  obtaining  an 
optimal  IBC  model.  It  was  also  shown  that  the  resulting  IBC  model  could  be  readily 
incorporated  into  an  existing  CEM  code  to  assess  the  performance  of  the  corrugated 
coating  when  mounted  on  complex  platforms. 

C.  FOLLOW-UP  STATEMENT: 

It  is  well  known  that  antenna  characteristics  such  as  input  impedance,  mutual 
coupling  and  radiation  pattern  can  be  drastically  altered  by  the  platform  that  supports  the 
antenna  structure.  In  order  that  platform  effects  are  properly  accounted  for,  the  platform 
must  be  considered  together  with  the  antenna  structure  during  simulation  and 
optimization.  This  usually  makes  the  design  process  a  laborious  and  time-consuming 


task.  It  is  therefore  desirable  to  develop  a  systematic  framework  where  such  design 
processes  are  highly  streamlined.  It  is  with  this  goal  that  we  have  carried  out  our 
research  under  the  H-infinity  Program.  To  accelerate  repeated  CEM  calculations  over 
large  parameter  spaces  such  as  frequency,  position  and  aspect  angle,  we  have  developed  a 
number  of  model-based  extrapolation/interpolation  algorithms  to  populate  the  antenna 
characteristics  data  based  on  a  very  sparse  set  of  CEM  calculations.  We  have  also 
developed  algorithms  to  extract  sparse  models  of  the  relevant  physics  embedded  in  the 
CEM  data  to  facilitate  design  and  synthesis. 

In  the  follow-on  research  program,  we  shall  move  on  to  explore  the  use  of  genetic 
algorithms  (GA)  for  shipboard  antenna  design,  placement  and  optimization.  Three 
problems  of  relevance  to  shipboard  antenna  design  will  be  addressed.  First,  we  propose 
to  investigate  the  optimal  design  of  electrically  small  wire  antennas  in  the  low  HF 
frequency  range,  where  the  mounting  platform  plays  a  prominent  role.  Second,  we 
propose  to  study  the  antenna  array  synthesis  and  placement  problem  in  the  presence  of 
mutual  coupling  and  platform  effects.  Third,  we  propose  to  investigate  the  shape  design 
of  microstrip  antenna  elements  for  broadband,  wide-angle  scanning  array  applications. 
While  GA  is  an  efficient  global  optimizer  ideally  suited  for  exploring  very  complex 
design  problems  such  as  those  outlined  above,  it  typically  requires  a  very  large  number  of 
cost  function  evaluations.  We  will  leverage  off  our  existing  research  to  minimize  the 
number  of  CEM  calculations  needed  for  cost  function  computation.  Particular  attention 
will  be  paid  to  the  close  coupling  between  GA  and  the  CEM  simulator  to  establish  a 
framework  in  which  the  design  process  can  be  streamlined.  The  potential  impact  of  the 
proposed  research  is  that  novel  design  concepts  will  be  uncovered  that  can  significantly 
outperform  conventional  designs  for  shipboard  antenna  problems. 

D.  REFERENCES: 

1.  R.  Roy,  A.  Paulraj  and  T.  Kailath,  “ESPRIT  -  a  subspace  rotation  approach  to 
estimation  of  parameters  of  cisoids  in  noise,”  IEEE  Trans.  Acoust.,  Speech,  Signal 
Processing,  vol.  ASSP-34,  pp.  1340-1342,  Oct.  1986. 

2.  L.  C.  Potter,  D.  Chiang,  R.  Carriere  and  M.  J.  Gerry,  “A  GTD-based  parametric 
model  for  radar  scattering,”  IEEE  Trans.  Antennas  Propagat.,  vol.  43,  pp.  1058- 
1067,  Oct.  1995. 


3.  J.  Song,  C.  C.  Lu  and  W.  C.  Chew,  “Multilevel  fast  multipole  algorithm  for 
electromagnetic  scattering  by  large  complex  objects,”  IEEE  Trans.  Antennas 
Propagat.,  vol.  45,  no.  10,  pp.  1488-1493,  Oct.  1997. 

4.  H.  T.  G.  Wang,  M.  L.  Sanders  and  A.  Woo,  “Radar  cross  section  measurement  data 
of  the  VFY  218  configuration,”  Tech.  Kept.  NAWCWPNS  TM-7621,  Naval  Air 
Warfare  Center,  China  Lake,  CA,  Jan.  1994. 

5.  R.  S.  Adve,  T.  K.  Sarkar,  D.  M.  Rao,  E.  K.  Miller,  and  D.  R.  Pflug,  “Application  of 
the  Cauchy  method  for  extrapolating/interpolating  narrow  band  system  responses,” 
IEEE  Trans.  Microwave  Theory  Tech.,  vol.  45,  pp.  837-845,  May  1997. 

6.  J.  Tsao  and  B.D.  Steinberg,  “Reduction  of  sidelobe  and  speckle  artifacts  in 
microwave  imaging:  the  CLEAN  technique,”  IEEE  Trans.  Antennas  Propagat.,  vol. 
AP-36,  pp.  543-556,  Apr.  1988. 

7.  T.  H.  Chu,  “Microwave  diversity  imaging  using  ASAR  approach,”  International 
IEEE  Antennas  and  Propagation  Symposium,  Salt  Lake  City,  UT,  July  2000. 

8.  S.  G.  Mallat  and  Z.  Zhang,  “Matching  pursuit  with  time-frequency  dictionaries,” 
IEEE  Trans.  Signal  Processing,  vol.  41,  pp.  3397-3415,  Dec.  1993. 

9.  A.  Kavak,  Vector  Propagation  Channel  Studies  for  Smart  Antenna  Wireless 
Communication  Systems,  Ph.D.  Dissertation,  The  University  of  Texas  at  Austin, 
Jan.  2000. 

10.  S.  M.  Rao,  D.  R.  Wilton,  and  A.  W.  Glisson,  “Electromagnetic  scattering  by 
surfaces  of  arbitrary  shape,”  IEEE  Trans.  Antennas  Propagat.,  vol.  30,  pp.  409-418, 
May  1982. 

11.  H.  Deng  and  H.  Ling,  “Fast  solution  of  electromagnetic  integral  equations  using 
adaptive  wavelet  packet  transform,”  IEEE  Trans.  Antennas  Propagat.,  vol.  AP-47, 
pp.  674-682,  Apr.  1999. 

12.  H.  Deng  and  H.  Ling,  “On  a  class  of  pre-defined  wavelet  packet  bases  for  efficient 
representation  of  electromagnetic  integral  equations,”  IEEE  Trans.  Antennas 
Propagat.,  vol.  AP-47,  pp.  1772-1779,  Dec.  1999. 

13.  Y.  Rahmat-Samii  and  E.  Michielssen,  Electromagnetic  Optimization  by  Genetic 
Algorithms,  John  Wiley  &  Sons:  New  York,  1999. 


E.  PUBLICATIONS: 

I.  LIST  OF  JOURNAL  ARTICLES  (ONR  supported  in  whole  or  in  part) 


14.  R.  Bhalla  and  H.  Ling,  “Near-field  signature  prediction  using  far-field  scattering 
centers  extracted  from  the  shooting  and  bouncing  ray  technique,”  IEEE  Trans. 
Antennas  Propagat.,  vol.  AP-48,  pp.  337-338,  February  2000. 

15.  Y.  Wang  and  H.  Ling,  “Adaptive  ISAR  image  construction  from  unevenly 
undersampled  data,”  IEEE  Trans.  Antennas  Propagat.,  vol.  AP-48,  pp.  329-331, 
February  2000. 

16.  T.  Su  and  H.  Ling,  “Determining  the  equivalent  impedance  boundary  condition  for 
corrugated  coatings  based  on  the  genetic  algorithm,”  IEEE  Trans.  Antennas 
Propagat.,  vol.  AP-48,  pp.  374-382,  March  2000. 

17.  C.  Ozdemir,  R.  Bhalla  and  H.  Ling,  “A  radiation  center  representation  of  antenna 
radiation  pattern  on  a  complex  platform,”  IEEE  Trans.  Antennas  Propagat.,  vol. 
AP-48,  pp.  992-1000,  June  2000. 

18.  T.  Su,  C.  Ozdemir  and  H.  Ling,  “On  extracting  the  radiation  center  representation 
of  antenna  radiation  patterns  on  a  complex  platform,”  Microwave  Optical  Tech. 
Lett.,  vol.  26,  pp.  4-7,  July  2000. 

1 9.  Y.  Wang  and  H.  Ling,  “A  frequency-aspect  extrapolation  algorithm  for  ISAR  image 
simulation  based  on  two-dimensional  ESPRIT,”  IEEE  Trans.  Geo.  Science  and 
Remote  Sensing,  vol.  38,  pp.  1743-1748,  July  2000. 

20.  T.  Su,  K.  Dandekar  and  H.  Ling,  “Simulation  of  mutual  coupling  effect  in  circular 
arrays  for  direction  finding  applications,”  Microwave  Optical  Tech.  Lett.,  vol.  26, 
pp.  331-336,  September  2000. 

21.  B.  Jiang,  T.  Su  and  H.  Ling,  “Frequency  interpolation  of  electromagnetic  scattering 
data  using  a  hybrid  model,”  Microwave  Optical  Tech.  Lett.,  vol.  27,  pp.  307-312, 
December  2000. 

22.  K.  R.  Dandekar,  H.  Ling  and  G.  Xu,  “The  effect  of  mutual  coupling  on  direction 
finding  in  smart  antenna  applications,”  to  appear  in  Elect.  Lett.,  December  2000. 

23.  T.  Su  and  H.  Ling,  “On  modeling  mutual  coupling  in  antenna  arrays  using  the 
coupling  matrix,”  to  appear  in  Microwave  Optical  Tech.  Lett.,  February  2001. 

24.  H.  Choo,  A.  Hutani,  L.  C.  Trintinalia  and  H.  Ling,  “Shape  optimization  of 
broadband  microstrip  antennas  using  the  genetic  algorithm,”  accepted  for 
publication  in  Elect.  Lett.,  November  2000. 

25.  H.  Deng  and  H.  Ling,  “A  wavelet-based  preconditioner  for  three-dimensional 
electromagnetic  integral  equations,”  accepted  for  publication  in  Elect.  Lett., 
November  2000. 


26.  Y.  Wang  and  H.  Ling,  “Efficient  radar  signature  prediction  using  a  frequency- 
aspect  interpolation  technique  based  on  adaptive  feature  extraction,”  submitted  for 
publication  in  IEEE  Trans.  Antennas  Propagat.,  June  1999. 

27.  T.  Su,  Y.  Wang  and  H.  Ling,  “A  frequency  extrapolation  technique  for  computing 
antenna-platform  radiation  problems,”  submitted  for  publication  in  IEEE  Trans. 
Antennas  Propagat.,  September  1999. 

28.  H.  Deng  and  H.  Ling,  “An  efficient  preconditioner  for  electromagnetic  integral 
equations  using  pre-defined  wavelet  packet  basis,”  submitted  for  publication  in 
IEEE  Trans.  Antennas  Propagat.,  November  1999. 

29.  K.  R.  Dandekar,  H.  Ling  and  G.  Xu,  “Experimental  study  of  mutual  coupling 
compensation  in  smart  antenna  applications,”  submitted  for  publication  in  IEEE 
Trans.  Communications,  August  2000. 

30.  L.  C.  Trintinalia  and  H.  Ling,  “Enhanced  triangular  patch  basis  functions  for 
electromagnetic  scattering  analysis,”  submitted  for  publication  in  IEEE  Trans. 
Antennas  Propagat.,  October  2000. 


H.  LIST  OF  CONFERENCE  PROCEEDINGS  (ONR  supported  in  whole  or  in  part) 

31.  H.  Choo,  T.  Su,  L.  C.  Trintinalia  and  H.  Ling,  “Shape  optimization  of  printed  patch 
structures  using  the  genetic  algorithm,”  URSI  National  Radio  Science  Meeting,  p. 
78,  Salt  Lake  City,  UT,  July  2000. 

32.  L.  C.  Trintinalia  and  H.  Ling,  “An  improved  triangular  patch  basis  for  the  method 
of  moments,”  International  IEEE  AP-S  Symposium,  pp.  2306-2309,  Salt  Lake  City, 
UT,  July  2000. 

33.  T.  Su  and  H.  Ling,  “A  model  for  the  active  element  pattern  of  array  elements 
including  both  mutual  coupling  and  platform  effects,”  International  IEEE  AP-S 
Symposium,  pp.  2130-2133,  Salt  Lake  City,  UT,  July  2000. 

34.  H.  Deng  and  H.  Ling,  “A  wavelet-packet  based  preconditioner  for  iterative  solution 
of  electromagnetic  integral  equations,”  International  IEEE  AP-S  Symposium,  pp. 
1810-1813,  Salt  Lake  City,  UT,  July  2000. 


ffl.  LIST  OF  RELATED  PRESENTATIONS 

35.  “Scattering  from  periodic  surfaces  -  analysis,  design  and  optimization,”  Lockheed 
Fort  Worth  Company,  Fort  Worth,  Texas,  April  10,  2000. 


36.  “Fast  postprocessing  algorithms  for  fast  CEM  solvers,”  2000  Electromagnetics  Code 
Consortium  Annual  Meeting,  Boeing  Phantomworks,  St.  Louis,  Missouri,  May  9, 
2000. 

37.  “Application  of  model-based  signal  processing  for  antenna  design,  placement  and 
optimization,”  Navy  H-Infinity  Program  Workshop,  US  Naval  Academy,  Annapolis, 
Maryland,  October  26, 2000. 


IV.  LIST  OF  THESES  AND  DISSERTATIONS 
Ph,D. 

I.  Y.  Kelly,  “The  multipath  fingerprint  method  for  wireless  E-91 1  location  finding,” 
May  2000. 

H.  Deng,  “Applications  of  wavelet  packet  bases  to  computational  electromagnetics 
and  radar  imaging,”  August  2000. 

M.S. 

J.  A.  Dobbins,  “Folded  conical  helix  antenna,”  May  2000. 

B.  Jiang,  “Model-based  interpolation  algorithms  for  electromagnetic  scattering  by 
complex  objects,”  May  2000. 

H.  Choo,  “Shape  optimization  of  printed  patch  structures  using  the  genetic 
algorithm,”  August  2000. 


V.  CONTRACTS  AND  GRANTS 

H.  Ling,  "Application  of  model-based  signal  processing  methods  to  computational 
electromagnetics  simulators,"  Office  of  Naval  Research,  December  1,  1997  - 
November  30, 2000. 

H.  Ling,  "MURI  center  for  computational  electromagnetics  research,"  Air  Force 
Office  of  Scientific  Research  (via  Univ.  of  Illinois),  December  15,  1995  -  December 
14, 2000. 

H.  Ling,  "Radar  image  enhancement,  feature  extraction  and  motion  compensation 
using  joint  time-frequency  techniques,"  Office  of  Naval  Research,  April  15,  1998  - 
September  30, 2001. 

H.  Ling,  "Electromagnetic  scattering  from  periodic  surfaces,"  Lockheed  Martin 
Corporate  Grant,  November  15,  1998  -  December  31, 1999. 

D.  T.  Jaffe  and  H.  Ling,  “High  index  grisms  for  mid-infrared  spectroscopy,” 
NASA,  June  1, 1999 -May  31, 2001. 


G.  Xu,  H.  Ling  and  H.  D.  Foltz,  “Development  of  wideband  vector  channel  models 
and  testbed  for  3rd  generation  wireless  mobile  systems,”  Texas  Advanced 
Technology  Program,  January  1, 2000  -  December  31, 2001. 


F.  INTERACTIONS/COLLABORATIONS  WITH  NAVY  SCIENTISTS: 

The  third  H-infinity  workshop  was  held  at  the  Naval  Academy  on  October  26,  2000. 
In  addition  to  reporting  our  progress  and  exchanging  ideas  with  the  other  team  members 
of  the  H-infinity  program,  new  collaborative  opportunities  were  explored.  We  continue 
to  interact  closely  with  Dr.  Victor  Chen  of  Naval  Research  Lab  on  a  separate  ONR 
program  in  applying  joint  time-frequency  processing  to  inverse  synthetic  aperture  radar 
imagery.  We  expect  to  generate  cross-fertilization  of  ideas  with  the  H-infinity  program 
since  a  good  physical  understanding  of  electromagnetic  scattering  and  radiation 
phenomena  is  the  basis  of  our  work  in  both  programs. 

G.  NEW  DISCOVERIES,  INVENTIONS,  OR  PATENT  DISCLOSURES: 

None. 

H.  HONORS  AND  AWARDS: 

None. 


APPENDIX 


Publications  Supported  by  ONR  Grant 


IEEE  TRANSACTIONS  ON  ANTENNAS  AND  PROPAGATION,  VOL.  48,  NO.  2,  FEBRUARY  2000 


Paper [14] 


Near-Field  Signature  Prediction  Using  Far-Field 
Scattering  Centers  Extracted  from  the  Shooting  and 
Bouncing  Ray  Technique 

Rajan  Bhalla  and  Hao  Ling 


Abstract — We  present  a  technique  to  predict  the  near-field  radar  cross 
section  (RCS)  of  a  target  using  the  far-field  scattering  centers  extracted 
from  the  shooting  and  bouncing  ray  (SBR)  technique.  The  results  gener¬ 
ated  using  this  methodology  are  verified  against  the  brute-force  SBR  cal¬ 
culations  for  near-field  scenarios.  It  is  demonstrated  that  this  technique  is  a 
fairly  accurate  and  very  efficient  way  to  generate  near-field  RCS  data  pro¬ 
vided  that  the  transceiver  is  not  very  close  to  the  target 

Index  Terms — Radar,  ray  shooting,  scattering  centers,  target  identiifica- 
tion. 


I.  Introduction 

For  many  radar  applications  the  transmitting  and  receiving  antennas 
are  located  in  the  near-field  zone  of  the  target.  Under  this  situation, 
the  distance  between  the  transmitting  antenna  and  the  target  may  not 
be  large  enough  to  treat  the  incident  field  arriving  at  the  target  as  a 
plane  wave.  Similarly,  the  distance  between  the  receiving  antenna  and 
the  target  may  not  be  large  enough  to  treat  the  scattered  field  arriving 
at  the  receiver  as  a  plane  wave.  Hence,  the  far-field  assumption  is  not 
valid  and  a  near-field  analysis  is  necessary.  One  of  the  existing  method¬ 
ologies  to  compute  near-field  radar  cross  section  (RCS)  is  to  use  the 
shooting  and  bouncing  ray  (SBR)  technique  [1],  [2].  In  the  SBR  tech¬ 
nique,  rays  are  launched  at  the  target  from  the  phase  center  of  the  trans¬ 
mitting  antenna  and  are  traced  according  to  the  law  of  geometrical  op¬ 
tics.  At  the  exit  point  of  each  ray,  a  ray-tube  integration  is  performed 
to  find  its  contribution  to  the  total  scattered  field  at  the  receiving  an¬ 
tenna,  taking  into  account  of  the  distance  from  the  ray  tube  to  the  re¬ 
ceiver  location.  Although  the  SBR  methodology  is  straightforward,  for 
typical  radar  applications  the  problem  involved  often  contains  a  large 
number  of  parameter  combinations.  For  example,  in  typical  missile  en¬ 
gagement  simulations,  the  scattered  field  is  of  interest  at  different  trans¬ 
ceiver  locations  and  for  different  antenna  patterns.  Each  combination  is 
a  new  problem  and  requires  that  the  whole  computation  be  carried  out 
each  time.  This  can  lead  to  an  inordinately  large  amount  of  computation 
time.  In  this  letter,  we  present  a  technique  to  predict  the  near-field  RCS 
of  a  target  by  taking  advantage  of  the  far-field  scattering  centers  ex¬ 
tracted  using  the  existing  capability  of  the  SBR  technique  [3],  [4].  The 
basic  idea  is  to  extract  the  far-field  three-dimensional  (3-D)  scattering 
centers  once  and  then  use  the  strengths  and  positions  of  the  extracted 
scattering  centers  to  reconstruct  the  near-field  RCS  at  different  trans¬ 
ceiver  locations.  Since  the  extraction  is  done  only  at  the  far-field  range 
and  the  reconstruction  can  be  carried  out  very  quickly,  this  method¬ 
ology  is  potentially  a  very  efficient  way  to  compute  near-field  RCS. 

II.  Near-Field  RCS  Prediction  Using  Far-Field 
Scattering  Centers 

It  is  well  known  that  the  scattered  far  field  from  an  electrically  large 
target  can  be  modeled  by  a  discrete  set  of  scattering  centers  on  the 

Manuscript  received  December  31,  1997. 

R.  Bhalla  was  with  the  University  of  Texas  at  Austin,  Austin,  TX  78712-1084 
USA.  He  is  now  with  SAIC-DEMACO,  Boston,  MA  01880  USA. 

H.  Ling  is  with  the  Department  of  Electrical  and  Computer  Engineering,  the 
University  of  Texas  at  Austin,  Austin,  TX  78712-1084  USA. 

Publisher  Item  Identifier  S  00 18-926X(00)0 1282-5. 


Target  Geometry 

Scattering  Center  Model 

•  •  Is 

a  r*.*i*£ 

•  v... 

* — > 

(a) 

Tx 


Fig.  1.  (a)  Three-dimensional  far-field  scattering  center  model  of  a  target,  (b) 

Near-field  RCS  prediction  based  on  far-field  scattering  centers. 

target  [see  Fig.  1(a)].  The  monostatic  3-D  scattering  centers  at  a  given 
incident  direction  can  be  sparsely  represented  by  a  set  of  numbers 
{A„,  R„},  where  An  is  the  strength  of  a  scattering  center  and  Rn  is 
its  3-D  location.  The  3-D  scattering  center  extraction  methodology  of 
a  target  using  the  SBR  technique  has  been  developed  recently  in  [3] 
and  [4],  To  predict  near- field  data  using  the  far-field  scattering  centers, 
we  argue  that  each  scattering  center  models  the  interaction  between 
the  incident  electromagnetic  wave  and  a  localized  portion  of  the 
target.  The  portion  of  the  target  that  gives  rise  to  each  scattering  center 
is  much  smaller  than  the  whole  target.  Consequently,  the  scattering 
phenomenon  represented  by  each  scattering  center  is  expected  to  be 
valid  at  a  distance  much  closer  than  the  usual  far-field  range  criterion 
of  L  =  2D2  /A,  where  D  is  the  maximum  target  dimension  and  A 
is  the  wavelength.  At  very  closed-in  ranges  the  scattering  centers 
themselves  will  be  defocused  and  this  interpretation  may  no  longer  be 
valid. 

Next,  we  provide  the  formula  to  reconstruct  the  monostatic  near- 
field  from  the  monostatic  far-field  scattering  centers.  To  reconstruct 
the  scattered  field,  Es  at  a  frequency  /  for  a  transceiver  located  at  Rt , 
we  sum  up  the  contribution  of  each  scattering  center  as  [see  Fig.  1(b)] 

E°(f,  Rt)  =  Y,  A"  -  Rn |) 

Notice  that  in  addition  to  the  phase  term  to  take  into  account  of  the 
near-field  effect,  we  also  use  a  range  correction  for  the  amplitude  of 
each  scattering  center.  The  scattered  field  defined  in  (1)  has  been  nor¬ 
malized  so  that  it  will  reduce  to  the  far-field  RCS  as  Rt  approaches 
the  far-field  range  [5],  [6].  If  the  antenna  pattern  is  not  isotropic,  the 
effect  of  the  antenna  pattern  can  also  be  accounted  for  in  (1).  In  ad¬ 
dition,  although  only  the  monostatic  scattering  is  considered  in  this 
paper,  the  methodology  is  extendible  to  the  more  general  bistatic  sce¬ 
nario.  It  should  also  be  noted  that  the  reconstruction  in  (1)  is  based  on 
single-bounce  scattering  mechanisms  and  is  not  very  accurate  for  scat¬ 
tering  centers  whose  contribution  comes  from  high  multibounce  mech¬ 
anism. 

III.  Results 

*  Two  numerical  examples  are  presented  to  demonstrate  the  accuracy 
of  this  technique.  The  reconstructed  near-field  RCS  from  far-field  scat¬ 
tering  centers  is  compared  to  the  brute-force  SBR  results  as  a  function 
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Fig.  2.  Comparison  of  the  near-field  RCS  computed  by  brute-force  SBR  and 
by  scattering  center  prediction  for  a  conical  target,  (a)  Target  geometry,  (b) 
Near-field  RCS  versus  range. 


of  range  of  the  transceiver.  The  antenna  pattern  of  the  transceiver  is  as¬ 
sumed  to  be  isotropic  and  horizontally  polarized.  The  first  example  is 
a  simple  canonical  target  comprising  of  two  triplates  mounted  on  an 
elongated  box  as  shown  in  Fig.  2(a).  The  monostatic  3-D  scattering 
centers  are  extracted  at  an  incident  azimuth  angle  of  — 15°  and  an  ele¬ 
vation  angle  of  15°.  The  far-field  criterion  for  this  3-m  target  is  760  m 
at  12.65  GHz.  The  near- field  RCS  is  computed  as  the  transceiver  is 
moved  closer  to  the  target  at  the  same  fixed  incident  angle  of  (- 15°, 
15°)  and  the  range  was  varied  from  320  to  15  m.  Fig.  2(b)  shows  the 
near-field  RCS  comparison  between  the  brute-force  computation  and 
the  scattering  center  reconstruction  plotted  on  a  logarithmic  scale  of 
the  range.  The  agreement  between  the  two  profiles  is  very  good,  even 
up  to  a  very  closed-in  range  of  50  m.  The  pattern  observed  in  the  RCS 
profile  is  consistent  with  the  beating  of  two  strong  scattering  centers  at 
the  two  triplate  positions. 

The  second  example  is  a  full-size  airplane.  The  monostatic  3-D  scat¬ 
tering  centers  are  extracted  at  an  incident  (AZ,  EL)  angle  of  (120°, 
—  1 8°).  The  far-field  distance  for  this  20-m  target  is  27  000  m  at  10  GHz. 
The  near-field  RCS  is  computed  as  the  transceiver  is  moved  closer  to 
the  target  at  the  same  fixed  incident  angle  of  (120°,  —  1 8°)  and  the  range 
was  varied  from  25  000  to  200  m.  Fig.  3  shows  the  near-field  RCS  com¬ 
parison  between  the  brute-force  SBR  computation  and  the  scattering 
center  reconstruction  plotted  on  a  logarithmic  scale  of  range.  Even 
though  this  target  has  some  multibounce  returns  from  inlet  ducts  whose 
reconstruction  will  not  be  very  accurate  using  scattering  center  predic¬ 
tions,  the  reconstructed  RCS  captures  the  salient  feature  of  the  brute- 


Fig.  3.  Comparison  of  the  near-field  RCS  computed  by  brute-force  SBR  and 
by  scattering  center  prediction  for  a  full-size  airplane. 


force  SBR  computation.  The  comparison  degrades  as  we  get  closer  to 
the  target.  In  terms  of  computation  time,  the  brute-force  computation 
takes  roughly  30  min/range  point  on  an  SGI  R4400  workstation.  The 
scattering  center  extraction  itself  can  be  carried  out  in  about  30  min. 
Once  the  scattering  centers  are  extracted,  however,  the  reconstruction 
can  be  done  in  seconds  at  all  the  range  values.  Hence,  the  scattering 
centers  offer  tremendous  computation  time  savings  for  the  computa¬ 
tion  of  near-field  RCS. 


IV.  Conclusion 

In  this  paper,  we  presented  an  efficient  technique  to  predict  near-field 
RCS  using  far-field  scattering  centers  extracted  from  the  SBR  tech¬ 
nique.  Such  a  scheme  allows  us  to  generate  the  near-field  RCS  in  real 
time  for  a  variety  of  input  parameters  such  as  range  and  antenna  pat¬ 
terns.  For  very  closed-in  ranges,  however,  this  scheme  is  expected  to 
break  down  as  the  scattering  centers  become  defocused.  In  addition, 
such  prediction  is  based  on  single-bounce  scattering  and  will  not  be  as 
accurate  for  multibounce  returns.  Although  only  near-field  monostatic 
scattering  is  considered  in  this  paper,  the  methodology  can  be  extended 
to  handle  more  general  bistatic  scenarios. 
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Adaptive  ISAR  Image  Construction  from  Nonuniformly 
Undersampled  Data 

Yuanxun  Wang  and  Hao  Ling 


Abstract — An  adaptive  approach  is  proposed  to  construct  ISAR  images 
from  nonuniformly  undersampled  data  in  the  angular  domain.  The  algo¬ 
rithm  uses  an  adaptive  scattering  feature  extraction  engine  in  place  of  the 
Fourier  transform  in  the  image  construction  procedure.  The  algorithm  en¬ 
tails  searching  and  extracting  out  individual  target  scattering  features  one 
at  a  time  in  an  iterative  fashion.  The  interference  between  different  target 
scattering  features  is  thus  avoided  and  a  clean  ISAR  image  without  the 
aliasing  effect  can  be  obtained.  The  algorithm  is  verified  by  constructing  the 
ISAR  image  from  the  chamber  measurement  data  of  the  model  VFY-218 
airplane. 

Index  Terms — Image  reconstruction,  synthetic  aperture  radar. 


I.  INTRODUCTION 

Constructing  an  inverse  synthetic  aperture  radar  (ISAR)  image  of  a 
target  requires  data  collection  in  both  the  frequency  and  angular  dimen¬ 
sions.  If  the  data  are  uniformly  sampled  and  the  sampling  rate  is  dense 
enough,  an  ISAR  image  can  be  obtained  by  using  a  two-dimensional 
(2-D)  fast  Fourier  transform  (FFT)  algorithm  [1].  In  this  paper,  we  ad¬ 
dress  the  case  when  the  angular  data  are  nonuniformly  undersampled. 
This  scenario  may  arise  in  real-world  data  collection  when  the  target 
is  fast  maneuvering  with  respect  to  the  radar  pulse  repetition  interval 
so  that  the  angular  look  on  the  target  by  the  radar  is  not  dense  enough 
to  satisfy  the  Nyquist  sampling  rate.  We  propose  an  algorithm  to  over¬ 
come  the  aliasing  effect  in  the  cross-range  dimension  and  construct 
ISAR  images  from  seriously  undersampled  data.  The  algorithm  uses 
an  adaptive  scattering  feature  extraction  engine  in  place  of  the  Fourier 
transform  in  the  image  construction  process.  The  original  concept  of 
adaptive  feature  extraction  was  proposed  in  [2]  and  [3].  It  has  been  ap¬ 
plied  to  ISAR  image  processing  in  the  joint  time-frequency  space  for 
resonant  scattering  mechanism  extraction  [4],  target  motion  compensa¬ 
tion  [5],  and  Doppler  interference  removal  [6].  In  contrast  to  the  Fourier 
transform,  where  the  signal  is  projected  to  all  the  image-domain  bases 
simultaneously,  the  adaptive  algorithm  searches  and  extracts  the  indi¬ 
vidual  target  scattering  features  one  at  a  time  in  an  iterative  fashion. 
When  applied  to  the  present  problem,  the  aliasing  error  caused  by  the 
interference  between  different  target  scattering  features  can  be  avoided. 
Therefore,  after  all  the  main  features  are  extracted,  they  can  be  dis¬ 
played  in  the  ISAR  image  plane  without  the  aliasing  effect.  We  verify 
this  algorithm  by  constructing  the  ISAR  image  using  the  chamber  mea¬ 
surement  data  of  the  model  VFY-218  airplane  [7].  It  is  found  that  a  rea¬ 
sonable  ISAR  image  can  be  constructed  from  seriously  undersampled 
data. 


U.  Adaptive  Feature  Extraction  (AFE)  Algorithm 

In  standard  ISAR  image  construction,  the  target  is  assumed  to  be  a 
collection  of  point  scattering  centers.  Under  the  small-angle  approx- 
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imation,  the  scattered  field  from  the  target  observed  as  a  function  of 
frequency  and  angle  can  be  written  as 


N 


E(f,9)  =  J2  0(xi,yi)e-2]kZi 
1=1 
N 


-  ^2  °iXi'Vi)e' 


2jkxi-2jkceyi 


i=l 


(1) 


where  0(xi,  pi)  is  the  amplitude  of  the  zth  scattering  center,  k  is  the 
free-space  wave  number,  and  kc  corresponds  to  the  wave  number  at 
the  center  frequency.  Xi  and  yl  denote  the  down  range  and  cross-range 
dimensions,  respectively.  We  assume  that  the  sampling  in  frequency 
is  uniformly  spaced  and  dense  enough  to  satisfy  the  Nyquist  criterion 
since  it  is  completely  controlled  by  the  radar.  Thus,  the  range  pro¬ 
file  versus  angle  data  can  be  generated  from  the  frequency- aspect  data 
by  applying  a  one-dimensional  (1-D)  Fourier  transform  along  the  fre¬ 
quency  dimension.  We  shall  denote  the  result  as  R(x,  8) 

N 

R(x,  d)  =  ^2  0{x"  yi^S*(x  ~  *i)e~2ikcVi<>-  (2) 

t=l 


In  the  above  expression,  Sx(x  —  Xi)  is  the  down-range  point  spread 
function  due  to  the  finite-length  frequency  domain  data.  Similarly,  the 
cross-range  information  can  also  be  obtained  from  angular  data  via  a 
1-D  Fourier  transform  of  R(x,  9)  along  the  angular  dimension.  The 
resulting  image  7(x,  y)  is 

I(x,y)  =  j  R(x,8)e2ik‘v°  d8 

N  r 

=  Y,0(xi,yi)SI(x-xi)  /  e2ik^-Vi)d8 
i=  1  ^ 

AT 

—  ^  ^  0{xi,  yi)Sx{x  —  Xi}Sy(y  -  y»)  (3a) 

i=i 

where 

Sy(y  -  Vi)  =  J  e2iM(w~Vi)  de  (3b) 

is  the  cross-range  point  spread  function  due  to  the  finite-length  angular 
domain  data.  If  the  data  are  sampled  densely  enough  such  that  the  nu¬ 
merical  integration  can  be  carried  out  accurately,  the  point-spread  func¬ 
tion  Sy  should  be  a  well-localized  function  with  its  peak  at  yi ,  while 
rapidly  decaying  away  from  the  peak.  The  resulting  image  I(x,  y)  will 
be  a  clean  image  with  good  indication  of  the  amplitudes  and  posi¬ 
tions  of  the  target  point  scattering  features.  However,  when  the  data 
are  undersampled,  the  numerical  integration  in  (3b)  will  result  in  large 
aliasing  error  that  shows  up  as  high  sidelobes  in  5y .  Consequently,  the 
constructed  image  will  contain  strong  interference  between  the  scat¬ 
tering  features.  This  effect  can  be  interpreted  as  the  loss  of  orthogo¬ 
nality  of  the  Fourier  bases  under  the  undersampled  condition. 

In  the  proposed  approach,  we  use  an  adaptive  feature  extraction  al¬ 
gorithm  in  place  of  the  Fourier  processing.  Instead  of  projecting  the 
signal  onto  all  the  exponential  bases  simultaneously,  we  search  out  the 
strongest  point  scattering  feature  in  the  cross-range  domain  and  remove 
it  from  the  original  signal.  Then  the  search  is  repeated  for  the  remainder 
signal  and  the  point-scattering  features  are  extracted  one  at  a  time  until 
the  energy  of  the  residue  signal  is  smaller  than  a  preset  threshold.  The 
search  procedure  is  carried  out  by  calculating  the  integral  in  (3b)  for  all 
points  in  cross  range  but  saving  only  the  maximum  value  and  position, 
i.e., 


[Bp,yp]  =  max  [/P(x,y)] 


(4) 


(dBsm) 


Fig.  1.  The  ISAR  image  constructed  at  30°  azimuth  from  the  original  data 
using  FFT. 


Range  (m) 


(dBsm) 


Fig.  2.  The  ISAR  image  constructed  at  30°  azimuth  from  randomly 
undersampled  data  using  Fourier  transform. 


where  p  denotes  the  stage  of  the  iterative  procedure.  The  remainder 
signal  is  produced  by  subtracting  out  the  pth  feature 

Rp+i(x,8)  =  Rp(x,9)  -  Bpe~2*kcVp6 .  (5) 

The  convergence  of  the  above  procedure  is  guaranteed  and  the  math¬ 
ematical  proof  is  given  in  [2].  The  advantage  of  such  an  iterative  pro¬ 
cedure  is  that  each  time  we  extract  out  the  strongest  feature,  we  also 
eliminate  its  interference  on  the  other  features.  It  should  be  noted  that 
nonuniform  sampling  is  a  prerequisite  to  ensure  that  there  is  no  ambi¬ 
guity  in  the  strongest  features  since  uniformly  undersampled  data  will 
result  in  multiple  positions  of  the  strongest  features.  For  simplicity,  the 
algorithm  is  repeated  for  each  range  cell  of  the  image.  A  2-D  algorithm 
in  frequency  and  aspect  can  also  be  developed,  if  the  search  is  carried 
out  for  both  the  range  and  cross-range  parameters.  After  all  the  features 
are  extracted  out,  we  can  construct  an  ISAR  image  using  the  amplitudes 
and  positions  of  the  point  scatterers. 

m.  Results  and  Discussion 

r  To  examine  the  applicability  of  the  algorithm  on  real  target  scattering 
data,  we  reconstruct  the  radar  image  of  a  model  (1:30  scale)  VFY-218 
airplane  using  undersampled  chamber  measurement  data  [7].  The  mea¬ 
surement  data  consist  of  an  aspect  window  from  10°  to  50°  and  a  fre¬ 
quency  range  from  8  to  16  GHz.  To  construct  an  ISAR  image,  we  first 
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Fig.  3.  The  ISAR  image  constructed  at  30°  azimuth  from  randomly 
undersampled  data  using  the  AFE  algorithm. 


polar  reformat  the  frequency-aspect  data  to  the  (Kx,I\y)  space.  The 
reformatted  data  consist  of  401  samples  in  Kx  and  438  samples  in  Ky . 
The  ISAR  image  is  first  generated  by  fast  Fourier  transform  (FFT)  and 
shown  with  the  airplane  overlay  in  Fig.  1.  The  point-scattering  features 
can  clearly  be  seen.  Next,  we  test  our  algorithm  by  generating  an  under¬ 
sampled  data  set  in  Ky .  This  is  approximately  the  same  as  undersam¬ 
pling  in  angle.  (Note  that  for  full  size  targets,  this  approximation  gets 
better.)  The  Nyquist  sampling  rate  requires  about  36  sampling  points 
in  Ky  and  we  randomly  select  only  24  out  of  the  438  points.  The  max¬ 
imum  sampling  interval  used  is  about  four  times  the  size  of  the  Nyquist 
sampling  interval.  Therefore,  serious  aliasing  will  occur  if  the  Fourier 
transform  algorithm  is  used,  as  shown  by  the  ISAR  image  displayed  in 
Fig.  2.  All  the  features  are  overlapped  with  sidelobe  noise  such  that  the 
point  scattering  features  on  the  airplane  can  no  longer  be  distinguished. 
Next,  we  apply  the  adaptive  feature  extraction  (AFE)  algorithm  to  each 
range  cell  of  the  image.  The  image  is  reconstructed  and  shown  in  Fig.  3. 
Comparing  Fig.  3  with  Fig.  1,  we  can  see  the  main  features  of  the  air¬ 
plane  in  Fig.  3  are  all  well  reconstructed  in  Fig.  1 .  We  do  observe  some 
noisy  spots  outside  the  target  at  the  lower  dynamic  ranges  in  Fig.  3. 
This  low-level  noise  occurs  at  about  25  dB  down  from  the  key  features 
and  presents  a  dynamic  range  limitation  of  the  present  AFE  algorithm. 
The  algorithm  has  also  been  tested  on  measured  data  from  in-flight  tar¬ 
gets  with  good  success. 
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Determining  the  Equivalent  Impedance  Boundary 
Condition  for  Corrugated  Coatings  Based  on  the 

Genetic  Algorithm 

Tao  Su,  Student  Member,  IEEE,  and  Hao  Ling 


Abstract — A  methodology  based  on  the  genetic  algorithm  (GA) 
is  proposed  to  determine  the  equivalent  impedance  boundary  con¬ 
dition  (IBC)  for  corrugated  material  coating  structures.  In  this  ap¬ 
proach,  rigorous  solutions  of  the  reflection  coefficients  at  a  number 
of  incident  angles  are  first  calculated  using  a  periodic  method  of 
moments  (MoM)  solver.  The  IBC  model  is  used  to  predict  the  re¬ 
flection  coefficients  at  the  same  observation  angles.  The  model  co¬ 
efficients  are  then  optimized  using  the  GA  so  that  the  difference 
between  the  approximated  and  the  MoM  predicted  reflection  coef¬ 
ficients  is  minimized.  The  GA  proves  efficient  in  obtaining  an  op¬ 
timal  IBC  model.  The  resulting  IBC  model  can  be  readily  incorpo¬ 
rated  into  an  existing  computational  electromagnetics  code  to  as¬ 
sess  the  performance  of  the  corrugated  coating  when  mounted  on 
complex  platforms. 

Index  Terms — Coatings,  genetic  algorithm,  gratings,  impedance 
boundary  conditions. 


I.  Introduction 

IT  is  well  known  that  the  impedance  boundary  condition 
(IBC)  approximation  is  an  efficient  way  to  model  complex 
structures  such  as  material  coatings  and  subskinline  features 
[l]-[3].  It  replaces  the  original  volumetric  structure  with  a 
surface  impedance  so  that  the  problem  dimension  is  reduced 
by  one.  Thus,  large  savings  in  computational  resources  can 
be  achieved  in  the  analysis  of  the  original  problem.  However, 
to  determine  a  simple  IBC  for  an  arbitrary  structure  that  is 
valid  over  a  wide  range  of  incident  angles,  polarizations  and 
frequencies  is  a  nontrivial  task.  In  this  paper,  we  set  out  to 
develop  a  methodology  to  determine  the  equivalent  IBC  model 
for  a  corrugated  coating  structure  backed  by  a  conducting 
surface  (see  Fig.  1).  The  corrugation  of  the  surface  is  assumed 
to  be  periodic  in  one  dimension  along  the  x-axis.  Of  interest  is 
an  IBC  model  that  is  valid  over  a  large  range  of  incident  angles 
in  both  the  0  and  ^  directions.  Our  objective  is  to  establish 
a  robust  methodology  such  that  the  resulting  IBC  model  can 
be  used  in  place  of  the  actual  coating  structure  in  subsequent 
analysis  and  design  involving  complex  platforms. 

The  problem  at  hand  is  difficult  since  the  scattering  charac¬ 
teristics  of  the  corrugated  surface  is  strongly  dependent  on  po- 
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Fig.  1 .  Scattering  from  a  corrugated  coating  structure  backed  by  a  conducting 
surface. 

larization  and  incident  angle.  The  standard  IBC  used  for  flat 
coatings  accounts  for  neither  the  aniostropic  nor  the  angular  be¬ 
havior  cf  the  scattering  characteristics  of  the  corrugated  surface. 

Some  improved  impedance  boundary  conditions  have  been  pro¬ 
posed  in  the  literature,  including  the  tensor  impedance  boundary 
condition  (TIBC)  [1]  and  the  generalized  impedance  boundary 
condition  (GIBC)  [2],  [4].  TIBC  usually  works  only  for  a  very 
limited  range  of  incident  angles.  GIBC  improves  the  accuracy 
of  the  IBC  model  by  including  higher  order  derivatives  of  the  f 

fields  on  the  surface.  However,  it  cannot  be  easily  implemented  f 

in  existing  MoM  solvers  since  it  requires  further  reformulation  f 

in  the  integral  equation.  A  resistive  boundary  condition  (RBC)  £ 

has  been  reported  that  works  well  over  large  incident  angles  for  * 

two-dimensional  (2-D)  planar  periodic  surfaces  [5].  However,  it 
is  limited  to  surfaces  with  very  small  periods.  Furthermore,  the 
choice  for  the  position  of  the  equivalent  impedance  surface  is 
not  obvious  for  the  corrugated  structure. 

Our  proposed  approach  to  this  problem  is  based  on  the  ge¬ 
netic  algorithm  (GA).  First,  we  compute  the  reflection  coeffi¬ 
cients  from  the  corrugated  surface  over  a  number  of  incident 
angles  and  polarizations  using  a  periodic  method  of  moments 
(MoM)  solver  [6].  The  resulting  reflection  coefficients  consti¬ 
tute  our  reference  data  base.  Next,  a  simple  periodic  IBC  model 
is  proposed  from  which  we  can  derive  an  expression  for  the  re¬ 
flection  coefficients.  In  the  GA  step,- we  optimize  the  model  co¬ 
efficients  so  that  the  difference  between  the  IBC-predicted  and 
the  MoM-predicted  reflection  coefficients  is  minimized.  GA 
searches  the  entire  parameter  space  in  a  way  similar  to  natural 
evolution  and  arrives,  after  many  generations,  at  the  best  param¬ 
eters  for  the  model. 
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(a)  (b) 

Fig.  2.  Equivalent  impedance  boundary  problem,  (a)  Original  structure,  (b) 
IBC  Approximation. 


This  paper  is  organized  as  follows.  In  Section  n,  the  MoM 
solution  of  the  problem,  the  IBC  model  formulation  and  the  GA 
optimization  are  discussed  as  the  steps  of  the  IBC  determina¬ 
tion.  Numerical  results  are  provided  in  Section  in  to  verify  the 
effectiveness  of  the  approach  in  a  number  of  corrugated  geome¬ 
tries.  A  3-D  scattering  example  is  also  given  to  demonstrate  the 
utility  of  the  resulting  IBC  model. 


n.  Methodology  for  Determining  the 
Optimal  IBC  Model 

In  this  section,  we  describe  the  proposed  methodology  for 
determining  the  equivalent  IBC  of  a  corrugated  coating  using 
the  GA.  In  the  first  step,  the  reflection  coefficients  from  the 
coating  are  computed  using  the  MoM  at  multiple  incident  angles 
to  serve  as  the  reference  data  of  the  model.  The  MoM  solution 
for  the  corrugated  coating  structure  in  Fig.  1  has  been  formu¬ 
lated  earlier  in  [6].  The  formulation  entails  dividing  one  cell  of 
the  grating  into  different  homogeneous  regions  according  to  the 
material  layers  as  shown  in  Fig.  2(a).  Boundary  integral  equa¬ 
tions  are  first  obtained  for  each  region.  Field  continuity  at  region 
interfaces  and  periodic  boundary  conditions  at  cell  boundaries 
are  then  enforced.  The  fields  in  the  top  half-space  are  expanded 
into  a  sum  of  Floquet  harmonics  and  are  matched  to  the  fields  in 
the  lower  region  so  that  the  reflection  coefficients  can  be  found. 

In  the  next  step,  a  periodic  IBC  model  is  proposed,  from 
which  we  can  derive  an  expression  for  the  reflection  coeffi¬ 
cients.  In  the  final  step,  the  optimal  parameters  for  the  IBC 
model  are  obtained  by  minimizing  the  mean  squared  error  be¬ 
tween  the  two  sets  of  reflection  coefficients  based  on  the  GA. 
These  steps  are  described  in  detail  below. 

A.  Periodic  IBC  Model 

The  equivalent  IBC  relating  the  tangential  electric  and  mag¬ 
netic  fields  for  a  planar  coated  surface  can  be  written  as  [1]: 


We  shall  adopt  this  model  for  the  corrugated  problem  due  to  its 
simplicity  and  usefulness  for  our  applications.  The  model  pa¬ 
rameters  will  then  be  optimized  to  emulate  the  properties  of  the 
exact  structure.  Note  that  since  the  corrugated  surface  exhibits 
anisotropic  scattering  characteristics,  the  equivalent  IBC  must 
also  in  general  be  anisotropic.  Therefore,  the  cross  impedance 
terms  Zxx  and  Zzz  are  kept  in  our  formulation  to  assess  their  im¬ 
portance.  The  boundary  impedance  Zzz,  Zzx,  Zxz,  and  Zxx  are 


Fig.  3.  Flow  chart  of  the  GA. 

in  general  functions  of  incident  angle  and  spatial  position.  For 
the  IBC  model  to  be  useful  for  subsequent  electromagnetic  anal¬ 
ysis,  however,  it  is  much  preferable  to  model  the  boundary  im¬ 
pedances  as  functions  of  spatial  position  only.  We  cannot  prove 
theoretically  the  existence  of  such  a  model  for  an  arbitrary  cor¬ 
rugated  structure.  Instead,  the  applicability  and  limitation  of  this 
approach  will  be  explored  numerically  in  Section  in. 

In  our  IBC  model,  the  periodic  grating  structure  with  period 
p,  as  shown  in  Fig.  2(a),  is  replaced  by  an  equivalent  impedance 
boundary  condition,  which  also  has  a  period  p,  as  illustrated  in 
Fig.  2(b).  Each  surface  impedance  term  can  be  expanded  into 
a  Fourier  series.  Since  the  cross  impedance  terms  Zzz  and  Zxx 
are  usually  very  small,  we  shall  treat  them  as  constants  and  only 
expand  the  impedances  Zzx  and  Zxz  as 

OO  oo 

Zzx  =  J2  ane~j3?nx>  Z*z=  J2  Ke~j^nx.  (2) 

n=— oo  n= — oo 

Therefore,  to  fully  describe  the  IBC  model,  we  must  determine 
the  Fourier  series  coefficients  {an}  and  {6n}. 

B.  Solution  to  the  Forward  Problem  of  Scattering  by 
the  IBC  Model 

Next,  we  derive  the  reflection  coefficients  resulting  from  the 
plane  wave  scattering  from  the  IBC  model  given  above.  Under 
plane  wave  incidence  where 

Ko  =  *o  sin  9  sin  <f>,  kxy0  =  -ko  cos  9  and 
Ko  =  *o  sin  6  cos  <f> 

each  component  of  the  tangential  electric  and  magnetic  fields  at 
the  impedance  surface  can  be  expanded  into  a  sum  of  Floquet 
harmonics  [7].  For  example,  the  tangential  electric  field  in  the  z 
direction  is  expanded  as 

oo 

Ez  =  Eie-jk’ax  +  .Ezne~jkmnX  O) 
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(a)  (b) 

Fig.  4.  Grating  geometry  for  (a)  triangular  groove  and  (b)  rectangular  groove. 

where 

I  i  27T 

n  —  ^xO  ”1“  ^ 

V 

is  the  propagation  constant  of  the  nth  order  harmonic 
along  the  x-direction.  This  Floquet  harmonic  is  a  reflected 
wave  with  propagation  constants  (kxn,kryn,kz),  where 
kyn  =  y/k%  -k\-  k*n  with  the  square  root  taken  as  positive 
real  or  negative  imaginary.  The  superscripts  i  and  r  denote 
respectively  the  incident  and  reflected  field  throughout  this 
paper.  The  harmonic  term  eju}te~^kzZ  is  suppressed  in  (3) 
and  y  is  set  to  zero  at  the  impedance  surface  for  convenience. 
Assuming  that  the  coefficients  of  the  higher  order  Floquet 
harmonics  are  negligible  and  applying  this  to  other  tangential 
field  components  we  get 

N 

F  =  Fie~jkx oX  -f  ^  F^e~jkxnX  (4) 

n~-N 

where  F  can  be  Ez,  Hz,  Ex  or  Hx ,  and  N  is  a  positive  integer. 
Substituting  (2)  and  (4)  into  (1)  and  matching  the  coefficients 
of  the  exponential  terms,  we  obtain  a  set  of  equations  relating 
Fzm  HXn  2nd  Hzn 

F\SnQ  +  Erzn 

N 

=  ZzzHTzn+  Y,  an-mHrxm+anH'x+ZzzHi6no 

m=  —  N 

(5) 

for  n  =  —N  to  N.  6n o  is  the  Kronecker  delta.  This  set  of  equa¬ 
tions  can  be  written  in  matrix  form  as 

El  uN  +  Erz  =  ZzzWz  +  AH;  +  Hi  a  +  ZzzH\  mn  (6) 

where 


Erz0  •• 

•  Kn)T 

H*  = 

[K-n  ■■■ 

HI o  • 

••  #r,iv]T 

h;  = 

ih:,-n  ••• 

Hrx o  • 

••  h^n}t 

ao 

•  •  •  a_7v  * 

*  a_2  N  1 

"a- at  ‘ 

A  = 

Q>N 

ao 

a~N  ,  a  = 

ao 

-  &2N 

aw 

ao  J 

-  aN  - 

Uyv  —  [u-AT  •**  Uq  •••  Un]T  ,  Un  =  SnQ. 


90  =  5  10  15  20  25  30  35  4045  50  55  60  65  70  75  80  85 
incident  angle 

(a) 


90 ~  5  10  15  20  25  30  35  40  45  50  55  60  65  70  75  80  85 
Incident  angle 

(b) 


90  -  5  10  15  20  25  30  35  4045  50  55  60  65  70  75  80  85 
Incident  angle 

(c) 

Fig.  5.  Comparison  of  reflection  coefficients  versus  angle  among:  1)  the 
exact  MOM  result;  2)  the  IBC  derived  from  GA;  and  3)  the  TTBC  result,  (a) 
H-poIarization.  (b)  V-polarization.  (c)  H-V  cross  polarization. 

r 

Following  the  same  steps,  we  also  arrive  at  the  relationship  for 

EXnt  HXnt  2nd  Uzn 

\  . 

bh;  =  e;  -  zxxwx  +  e{xun  -  zxxh'xun  -  Hib  (?) 


and 
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where 


-  b0 

b-N  *  *  * 

b-2N' 

“b-N " 

B  = 

bN 

bo 

b-N 

,  b  = 

bo 

-&2JV  **• 

bn 

bo  - 

-  bN  - 

For  a  plane  wave,  the  tangential  components  of  the  fields  in  the 
x-direction  Exy  Hx  can  be  expressed  in  terms  of  Ez  and  Hz  as 


where 
T  = 


Combine  (6)-(8),  the  matrix  relationship  between  the  incident 
and  reflected  fields  is  written  as 

-l 


.»;]  -  ([o 
■([ 


-Z„I 

B 

0  a 

ujv  —  Zzzun 


-  ’?  -ti][T1) 

[T*]  + 


-Uat  -ZxxUn  \ 

0  b  J ) 


El 

Hi 


(9) 


In  (9),  the  z  components  of  the  scattered  field  can  be  directly 
related  to  those  of  the  incident  field.  We  can  therefore  define 
the  reflection  coefficients  for  the  different  polarizations  as 

fiTM-TM  =  ^|  ]  RTM- 


\H'=0 


-TE  _  f?o  Hrzn 

~  E\ 


jfTE-TE  — 


Hi 


Ei=  0 


pTE-TM  _  Erxn 
^  “  rioHi 


H'=  0 


Ei=0 


(10) 

To  summarize,  the  reflection  coefficients  from  the  IBC  model 
can  be  calculated  by  using  (9)  and  (10)  if  an  and  bn  are  given. 
This  relationship  is  utilized  by  the  G A  to  calculate  the  reflection 
coefficients  for  a  given  sample  of  {an}  and  {bn}  and  compare 
them  with  the  reference  data  obtained  from  the  MoM  solution 
to  optimize  the  model  parameters. 

Two  comments  are  in  order.  First,  although  the  reflection  co¬ 
efficients  are  derived  for  the  TEZ/TMZ  polarizations,  they  can 
be  easily  transformed  to  the  more  conventional  vertical  (V)/hor- 
izontal  (H)  polarizations  with  respect  to  the  surface.  Second,  we 
have  assumed  the  position  of  the  IBC  surface  is  at  y  =  0  in  the 
above  formulation.  However,  if  the  IBC  surface  is  at  the  position 
y  =  -d  [as  shown  in  Fig.  2(b)],  a  factor  of  e~^2kyd  should  be 
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Fig.  6.  Effect  of  incorporating  the  cross  impedance  terms  in  the  IBC  model 
generated  from  GA,  (a)  H-polarization,  (b)  V-polarization. 


multiplied  to  each  of  the  reflection  coefficients  to  arrive  at  the 
correct  answer.  It  is  well  known  that  in  general  layered  coating 
problems,  there  is  not  a  preferred  position  for  the  impedance 
surface.  The  solution  can  sometimes  be  improved  by  applying 
the  IBC  at  a  position  other  than  a  natural  interface  in  the  struc¬ 
ture  [4].  Therefore,  by  including  the  position  of  the  IBC  surface 
as  an  additional  tuning  parameter  in  our  IBC  model,  we  can  fur¬ 
ther  improve  the  accuracy  of  the  model. 

C.  Genetic  Algorithm  to  Determine  the  Optimal 
IBC  Parameters 

In  the  GA,  the  parameters  to  be  optimized  are  first  encoded 
into  binary  form.  A  set  of  the  encoded  parameters  is  known  as  a 
chromosome.  The  basic  idea  of  GA  is  to  generate  a  pool  of  chro¬ 
mosomes,  discard  the  bad  ones,  keep  the  best  ones  and  let  them 
evolve  to  produce  better  chromosomes.  The  evaluation  of  each 
chromosome  is  performed  by  a  cost  function  which,  in  this  case, 
is  chosen  to  be  the  mean-squared  error  between  the  MoM  com¬ 
puted  reflection  coefficients  and  those  solved  using  (9)  and  (10) 
with  the  (Zzz,  an,bn,  ZXX)d )  parameters  decoded  from  the  cor¬ 
responding  chromosome.  Chromosomes  in  the  pool  are  ranked 
according  to  the  cost  function.  The  best  ones  are  selected  in  pairs 
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Incident  angle 


Incident  angle 


Fig.  7.  Effect  of  increasing  the  model  order  in  the  IBC  model  generated  from  GA.  (a)  Amplitude,  H-polarization.  (b)  Phase,  H-polarization.  (c)  Amplitude, 
V-polarization.  (d)  Phase,  V-poIarization. 


to  act  as  parents  of  the  next  generation.  Reproduction  of  children 
chromosomes  is  based  on  specific  rules  of  heredity  and  muta¬ 
tion.  The  process  of  selection  and  reproduction  is  repeated  until 
a  set  of  satisfactory  parameters  is  found  or  the  generation  limit 
is  reached.  The  flow  chart  of  the  GA  is  shown  in  Fig.  3.  Detailed 
discussion  of  the  GAs  can  be  found  in  [8]. 

In  the  IBC  model  for  the  corrugated  coating,  the  parameters 
to  be  optimized  are  the  coefficients  of  the  Fourier  expansion  an 
and  6n,  the  cross  impedances  Zzz,  Zxx  and  the  position  of  the 
impedance  surface  d.  Each  of  the  parameters  an,bn,Zzz  and 
Zxx  consists  of  a  real  part  and  an  imaginary  part.  We  assume  a 
symmetric  structure  so  that  an  =  a_n,6n  =  6_n.  For  an  ap¬ 
proximation  truncated  to  the  Nth  order,  the  total  number  of  real 
numbers  is  4N  4-  9.  The  number  of  bits  contained  in  each  pa¬ 
rameter  B  is  adjustable.  If  B  is  too  large,  the  convergence  of  GA 
will  be  slow.  If  B  is  too  small,  the  accuracy  of  the  calculation 
will  suffer.  In  the  examples  given  in  this  paper,  we  choose  B  =  8 
to  be  efficient  in  both  speed  and  accuracy.  In  order  to  encode  the 
unknown  parameters  into  binary  form,  the  minimum  and  max¬ 
imum  possible  values  of  each  parameter  are  required.  For  ex¬ 
ample,  the  values  for  the  real  and  imaginary  parts  of  the  Fourier 
series  an  and  bn  are  estimated  to  be  in  the  range  from  —3t}q 
to  3t7q >  which  is  found  to  be  reasonable  in  the  numerical  ex¬ 


amples.  Thus,  the  eight-bit  binary  00000000  denotes  — 3770  and 
11111111  represents  3t/o.  Zzz  and  Zxx  are  relatively  small  and 
their  real  and  imaginary  parts  vary  from  — O.I770  to  O.I770.  The 
distance  d  is  limited  between  the  upper  and  lower  boundaries  of 
the  coating  so  that  the  resulting  IBC  model  will  not  cause  any 
ambiguity  in  its  applications. 

In  the  beginning  of  the  GA,  a  number  of  chromosomes 
are  randomly  generated.  Each  chromosome  is  decoded  into 
parameters  Zzz,an,bn,  Zxx  and  d.  The  reflection  coefficients 
are  then  computed  using  (9)  and  (10).  The  cost  function  gives 
the  mean-squared  error  between  these  reflection  coefficients 
and  their  corresponding  MoM  solution 

Cost  =  52  |flPl-P2  -  Rol~P*\2 
o,4>  Pi  ,f2 

where  Rq  denotes  the  MoM  solution  of  the  reflection  co¬ 
efficients  at  a  specific  observation  angle  ( Q,4> )  and  P\,P2 
is  the  polarization  TEZ/TMZ  (or  V/H).  The  fitness  value 
of  each  chromosome  is  given  by 

f(Ci)  =  cx-  c2Cost(Ci) 
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Fig.  8.  Higher  order  reflection  coefficients  simulated  by  the  IBC  model,  (a)  Oth 


where  c\  and  c2  are  constants  and  C{  is  the  zth  chromosome 
in  the  population.  This  fitness  value  is  used  in  ranking  the 
chromosomes  and  selecting  of  parents  for  the  next  generation 
[8]— [1 1].  There  are  several  standard  ways  of  selection.  In 
this  paper,  the  roulette  wheel  selection  is  used  in  which  the 
probability  of  each  chromosome  to  be  selected  is  proportional 
to  its  fitness  value. 

After  two  chromosomes  are  selected,  they  mate  to  generate 
children.  This  is  realized  by  the  process  of  crossover  in  which  a 
break  point  is  randomly  chosen  in  the  chromosomes  and  the  two 
chromosomes  are  switched  at  that  point.  Mutation  is  imposed 
at  this  point  so  that  new  genes  appear  in  the  next  generation. 
The  mutation  rate,  which  is  the  portion  of  bits  to  be  randomly 
changed,  is  also  an  important  parameter  in  GA.  Experiments 
show  that  a  mutation  rate  of  5-8%  is  often  efficient  in  the  cal¬ 
culation. 

The  process  of  evaluation,  selection,  and  reproduction  is  re¬ 
peated  until  a  desired  mean-squared  error  is  achieved  or  a  max¬ 
imum  generation  is  reached.  For  a  population  of  400  chromo¬ 
somes,  the  Oth  order  IBC  (i.e.,  N  —  0)  takes  10-20  generations 
to  converge  to  the  optimum  while  the  second-order  IBC  takes 
200-400  generations. 
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order,  H-pol.  (b)  First  order,  H-pol.  (c)  Oth  order,  V-pol.  (d)  1st  order,  V-pol. 


III.  Numerical  Examples 

In  this  section,  some  examples  are  presented  to  demonstrate 
the  effectiveness  of  the  method.  The  first  example  is  a  deep  tri¬ 
angular  grooved  grating  with  relatively  small  period.  The  ge¬ 
ometry  of  one  cell  of  the  grating  is  shown  in  Fig.  4(a)  where 
the  period  p  =  0.067Ao  and  hi  =  0.22Ao,/&2  =  0.017Ao. 
The  coating  material  is  MagRAM  with  material  constants  er  = 
14.35  -  j0.28  and  fir  =  1.525  -  j 1.347  at  the  frequency  of  10 
GHz.  Seventeen  observation  angles  are  selected  which  include 
normal  incidence  and  the  combination  of  9  —  20° ,  40° ,  60° ,  80° 
and  (j)  =  0°,  30°,  60°,  90°.  In  this  example,  we  set  Zxx  and  Zzz 
in  (1)  to  zero  and  N  —  0  in  (4)  to  make  the  model  compa¬ 
rable  with  TIBC.  The  co-polarization  reflection  coefficients  for 
the  H-pol  and  V-pol  incidence  are  plotted  in  Fig.  5(a)  and  (b), 
respectively,  and  the  H-V  cross-polarization  reflection  coeffi¬ 
cients  are  plotted  in  Fig.  5(c).  In  the  figures,  the  z-axis  is  divided 
'into  sections  of  different  incident  angle  9 ,  which  varies  from  5° 
to  85°  in  steps  of  10°.  In  each  of  the  9  section,  the  grating  angle 
(j)  varies  from  0°  to  90°.  The  matching  of  the  reflection  coeffi¬ 
cients  between  the  GA  approach  and  MoM  solution  is  good  at 
most  incident  angles,  even  near  grazing  incidence.  The  value  of 
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Fig.  9.  Error  of  the  IBC  model  as  a  function  of  material  constant  e  =  e'  — je ". 


d  is  found  to  be  0.027Ao  from  the  tip  of  the  groove  or  0.21Ao 
above  the  ground  plane.  The  TIBC  result  is  also  generated  by 
using  the  reflection  coefficient  at  normal  incidence  to  derive  the 
equivalent  boundary  condition.  The  impedance  surface  is  placed 
at  the  plane  of  the  conductor  backing.  It  can  be  seen  that  the 
TIBC  results  deviate  significantly  from  the  reference  solution 
away  from  normal  incidence.  With  the  same  complexity  of  the 
boundary  condition,  GA  achieves  a  much  better  matching  be¬ 
cause  more  observation  points  are  used  in  the  modeling  and  be¬ 
cause  of  the  additional  degree  of  freedom  in  the  position  of  the 
IBC  surface. 

Next,  we  compare  the  IBC  approximations  with  and  without 
the  cross  impedance  terms  Zxx  and  Zzz.  The  structure  is  a  rect¬ 
angular  groove  as  shown  in  Fig.  4(b)  with  a  period  p  =  0.17Ao 
and  a  groove  depth  of  hi  =  hi  —  0.042A0.  The  material  con¬ 
stants  are  er  =  8.3  -  j2A  and  pr  -2-  j0.9  at  10  GHz.  The 
same  observation  points  are  used  as  in  the  previous  example. 
The  0th  order  (N  =  0)  IBC  is  determined  and  the  comparison 
between  the  reflection  coefficients  is  illustrated  in  Fig.  6.  With 
the  cross  impedance  included,  the  accuracy  of  the  approxima¬ 
tion  is  improved. 

In  the  third  example,  the  IBC  of  different  orders  are 
obtained  for  the  structure  shown  in  Fig.  4(b)  where  p  = 
0.42Ao,/ii  =  0.25Ao  and  =  0.17Ao.  The  coating  material 
is  the  same  as  that  in  example  2  but  the  period  is  much  larger 
and  the  groove  is  deeper.  The  reflection  coefficients  pre¬ 
dicted  by  the  Oth-order  and  second-order  model  are  plotted 
in  Fig.  7.  While  the  approximation  by  the  Oth-order  model  is 
fairly  satisfactory,  the  second-order  model  further  improves 
the  result  and  the  matching  is  better  at  most  incident  an¬ 
gles.  The  price  of  the  improvement  is  the  computation  time. 
For  the  Oth-order  modeling,  it  takes  only  a  few  minutes  for 
the  GA  to  converge  while  the  second-order  IBC  takes  more 
than  1  h  on  an  SGI  02  workstation  (RlOOOO/155  MHz). 
Another  consequence  as  the  order  of  the  model  is  increased 
is  that  the  resulting  IBC  will  show  more  spatial  variation. 
This  implies  that  when  the  IBC  model  is  utilized  in  subse¬ 
quent  analysis  using  numerical  electromagnetics  solvers,  the 
impedance  surface  must  be  divided  more  finely  to  faithfully 


Fig.  10.  Geometry  of  the  comer  reflector. 


describe  the  IBC.  This  will  lead  to  a  higher  computation  cost. 
Thus,  the  higher  order  model  is  not  recommended  unless  the 
Oth-order  one  is  intolerable  or  the  period  is  large  compared 
to  the  wavelength.  Generally  speaking,  the  IBC  model  can  be 
improved  by  increasing  the  model  order  whether  or  not  the 
cross-impedance  terms  exist.  But  with  the  cross-impedance 
terms,  the  required  model  order  is  usually  smaller  than  that 
without  them. 

We  further  investigate  a  structure  with  a  period  larger 
than  half  a  wavelength,  which  results  in  higher  order  Ro¬ 
quet-mode  reflection  at  some  incident  angles.  The  rectangular 
groove  shown  in  Fig.  4(b)  has  a  period  p  =  0.85Ao  and 
hi  —  h2  =  0.42Ao.  The  material  constant  are  er  —  10.5  —  j'2.2 
and  pr  =  2  -  j0.3.  The  0th  and  1st  order  reflection  coefficients 
are  plotted  in  Fig.  8(a)  and  (b),  respectively.  It  is  shown  that  the 
Roquet  modes  are  also  well  characterized. 

In  the  final  example,  we  investigate  the  limitation  of 
the  IBC  model.  We  consider  a  triangular  groove  shown  in 
Fig.  4(a)  with  p  =  0.33Ao  and  hi  =  /i2  =  O.O8A0.  The 
optimal  IBC  model  is  found  using  the  GA  for  different 
coating  materials.  A  second-order  IBC  model  with  cross 
impedance  terms  is  used  and  the  optimal  model  parameters 
are  determined  by  running  GA  to  convergence.  After  the 
model  is  found,  the  root  mean  squared  (RMS)  error  of 
the  IBC-predicted  reflection  coefficients  over  the  selected 
angles  are  computed.  Fig.  9  shows  the  RMS  error  of  the 
optimal  IBC  model  as  a  function  of  efr  and  e”>  which 
are  the  real  and  imaginary  parts  of  the  coating  relative 
permittivity  er .  pr  is  set  to  one  for  all  the  coatings.  We 
observe  that  the  IBC  model  works  best  for  high-contrast 
high-loss  materials.  For  low-contrast  or  low-loss  materials, 
the  model  error  can  be  large.  This  behavior  is  very  similar 
to  conventional  IBC  models  for  planar  coatings. 

We  now  apply  our  derived  IBC  model  to  a  three-dimen¬ 
sional  (3-D)  scattering  problem.  Consider  the  comer  reflector 
as  shown  in  Fig.  10.  The  monostatic  radar  cross  section 
(RCS)  is  calculated  for  both  the  uncoated  reflector  and  that 
coated  with  the  MagRAM  structure  described  in  example 
1.  The  groove  of  the  coating  is  either  parallel  or  perpen¬ 
dicular  to  the  incident  direction.  Both  cases  are  computed 
for  comparison.  Note  that  the  solution  for  such  a  structure 
is  very  complicated  if  we  try  to  use  the  exact  MoM  for¬ 
mulation.  Instead,  we  use  the  Oth-order  impedance  boundary 
condition  obtained  from  example  1  to  replace  the  corrugated 
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Fig.  11.  RCS  of  the  comer  reflector  (/  =  10  GHz),  (a)  H-polarization.  (b) 
V-polarization. 


absorber.  We  assume  that  the  size  of  the  plate  remain  the 
same  after  the  IBC  replacement.  The  RCS  is  computed  using 
FISC  [12],  which  is  a  3-D  MoM  code  based  on  the  fast 
multipole  method  [13].  Comparisons  of  the  RCS  at  several 
elevation  angles  are  shown  in  Fig.  1 1  for  both  the  H-  and 
V-polarizations.  The  result  shows  the  effect  of  coating,  which 
lowers  the  overall  RCS  level  for  both  polarizations.  We  fur¬ 
ther  observe  that  a  30-dB  RCS  reduction  can  be  achieved  for 
both  polarizations  over  the  range  of  elevation  angles  from 
20°  to  75°  if  the  grating  is  oriented  parallel  to  the  incident 
wave. 


IV.  Conclusion 

In  this  paper,  an  impedance  boundary  condition  model  is 
derived  based  on  the  GA  to  approximate  arbitrary  corrugated 


coating  structures  in  scattering  problems.  The  periodic  struc¬ 
ture  is  replaced  by  a  periodic  IBC  on  a  virtual  surface.  The 
boundary  impedance  and  the  position  of  the  surface  are  opti¬ 
mized  by  matching  the  reflection  coefficients  to  the  rigorous 
numerical  solution  at  a  number  of  incident  angles.  Similar 
to  traditional  IBC  models,  this  approach  is  most  effective 
when  the  coating  material  is  high  loss  and  of  high  con¬ 
trast.  The  resulting  IBC  model  generated  by  this  algorithm 
can  be  incorporated  into  an  existing  computational  electro¬ 
magnetics  code  to  assess  the  performance  of  the  corrugated 
coating  when  mounted  on  complex  platforms. 

Compared  with  other  IBC  approaches,  the  method  de¬ 
scribed  above  has  several  advantages.  First,  the  boundary 
impedance  is  assumed  to  be  anisotropic  so  that  the  same 
model  can  be  applied  to  oblique  incidence  from  any  arbitrary 
angles.  Second,  it  is  possible  to  build  in  spatial  variation  of 
the  boundary  impedance  by  adjusting  the  number  of  terms 
used  in  the  Fourier  series  expansion.  By  using  more  terms, 
the  IBC  model  can  be  made  more  accurate.  In  addition,  the 
position  of  the  impedance  surface  can  also  be  optimized. 
By  solving  for  the  best  position  of  the  impedance  surface 
as  one  of  the  model  parameters,  the  accuracy  of  the  model 
can  be  improved. 

Numerical  experiments  show  that  the  IBC  approximation 
can  be  improved  if  some  of  the  parameters  of  the  GA  are 
carefully  chosen.  These  parameters  include  the  incident  an¬ 
gles  at  which  the  rigorous  solution  is  obtained,  the  range 
of  each  model  parameter,  and  the  mutation  rate,  etc.  The 
GA  can  also  be  accelerated  with  carefully  chosen  parame¬ 
ters  and  a  well-designed  cost  function. 
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A  Radiation  Center  Representation  of  Antenna 
Radiation  Patterns  on  a  Complex  Platform 

Caner  Ozdemir,  Member,  IEEE,  Rajan  Bhalla,  and  Hao  Ling,  Fellow,  IEEE 


Abstract — A  sparse  model  of  the  antenna  radiation  pattern  on 
a  complex  platform  is  presented.  This  representation  is  based  on 
a  point  radiator  model  that  describes  the  radiation  pattern  by  a 
collection  of  radiation  centers  on  the  platform.  The  methodology 
for  obtaining  the  radiation  center  model  is  presented.  It  entails 
first  generating  the  three-dimensional  (3-D)  antenna  synthetic 
aperture  radar  (ASAR)  imagery  of  the  platform  and  then  param¬ 
eterizing  the  resulting  image  by  a  collection  of  point  radiators  via 
the  CLEAN  algorithm.  It  is  shown  that  once  such  a  representation 
is  obtained,  we  can  reconstruct  and  extrapolate  antenna  radiation 
patterns  over  frequencies  and  aspects  with  good  fidelity,  thus 
achieving  high  data  compression  ratio.  Furthermore,  it  is  shown 
that  the  resulting  radiation  center  information  can  be  used  to 
pinpoint  cause-and-effect  in  platform  scattering  and  provide 
important  guidelines  for  reducing  platform  effects. 

Index  Terms — Antenna  radiation  patterns,  scattering  centers. 


I.  Introduction 

IT  IS  well  known  that  the  platform  structure  that  supports  an 
antenna  can  dramatically  alter  its  radiation  characteristics. 
Recently,  we  developed  an  imaging  algorithm,  termed  antenna 
synthetic  aperture  radar  (ASAR)  imaging,  to  pinpoint  the 
locations  of  secondary  scattering  off  a  mounting  platform  from 
the  antenna  radiation  data  [1].  Our  approach  is  similar  to  the 
inverse  synthetic  aperture  radar  (ISAR)  concept.  Contrary  to 
conventional  ISAR  imaging,  a  key  complication  of  the  ASAR 
imaging  scenario  is  that  the  antenna  is  located  in  the  near 
field  of  the  platform.  It  was  shown  that  under  the  small-angle 
approximation  and  single-bounce  assumption,  a  Fourier  trans¬ 
form  relationship  exists  between  multifrequency,  multi-aspect 
radiation  data  and  the  three-dimensional  (3-D)  positions  and 
strengths  of  the  secondary  scatterers  on  the  platform.  Therefore, 
a  3-D  image  showing  the  spatial  locations  of  platform  scattering 
can  be  constructed  via  Fourier  inversion  of  the  radiation  data. 
This  concept  was  demonstrated  using  the  computed  radiation 
data  from  the  code  Apatch  [2],  which  employs  the  shooting 
and  bouncing  ray  (SBR)  technique  [3]-[6].  Fig.  1  shows  the 
CAD  model  of  an  airplane  with  a  small  vertical  dipole  antenna 
placed  above  the  cockpit.  The  entire  airplane  is  assumed  to 
be  perfectly  conducting.  By  collecting  the  multiple  frequency 
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(10-GHz  with  0.56-GHz  bandwidth)  and  multiple  aspect 
(3.33°  azimuth  window  and  1.72°  elevation  window  about  the 
nose)  radiation  pattern,  a  3-D  image  can  be  formed  from  the 
ASAR  algorithm.  Shown  in  Fig.  2(a)  is  the  3-D  ASAR  image 
viewed  along  various  horizontal  cuts  through  the  airplane. 
The  key  scattering  locations  near  the  nose  and  wings  of  the 
airplane  can  be  clearly  identified.  This  is  very  similar  to  an 
ISAR  image  in  which  the  strong  scattering  locations  on  the 
target  are  displayed  to  facilitate  signature  reduction  and  target 
identification  applications. 

In  this  paper,  we  extend  the  work  in  [1]  by  extracting  a  sparse 
representation  of  the  ASAR  imagery.  This  is  motivated  by  the 
observation  that  similar  to  an  ISAR  image,  ASAR  images  ex¬ 
hibit  strong  point-scatterer  like  behavior.  Since  it  is  well  known 
that  the  backscattered  signature  can  be  modeled  by  a  very  sparse 
set  of  scattering  centers  on  the  target  [6]-[9],  we  believe  it  is 
also  possible  to  parameterize  the  radiation  patterns  of  an  an¬ 
tenna  mounted  on  a  complex  platform  via  a  set  of  “ radiation 
centers”  Using  an  approach  similar  to  our  recent  work  on  scat¬ 
tering  center  extraction  [6],  we  utilize  the  CLEAN  algorithm 
[10],  [11]  to  carry  out  the  radiation  center  extraction  process 
from  an  ASAR  image.  It  is  shown  that  once  such  a  representa¬ 
tion  is  obtained,  we  can  reconstruct  and  extrapolate  antenna  ra¬ 
diation  patterns  over  frequencies  and  angles  with  good  fidelity, 
thus  achieving  high  data  compression  ratio.  Furthermore,  it  is 
shown  that  when  coupled  with  the  SBR  simulation  engine,  the 
resulting  radiation  center  information  can  be  used  to  pinpoint 
cause-and-effect  in  platform  scattering  and  provide  important 
guidelines  for  reducing  platform  effects.  :  .  v .  , 

This  paper  is  organized  as  follows.  In  Section  n,  we  first 
review  the  ASAR  image  formation  algorithm  and  present  the 
radiation  center  model  of  the  antenna  radiation  data.  Next,  we 
apply  the  CLEAN  algorithm  to  carry  out  the  extraction  process 
and  examine  the  sparseness  of  the  model  and  the  fidelity  of  the 
reconstructed  data  in  frequency  and  aspect.  In  Section  IH,  we 
demonstrate  two  utilities  of  the  radiation  center  representation. 
In  the  first  example,  we  show  that  the  antenna  radiation  pattern 
can  be  extrapolated  over  aspect  due  to  the  angular  stability  of 
the  radiation  center  model.  Consequently,  it  is  possible  to  rep¬ 
resent  the  complex  radiation  pattern  over  all  observable  angles 
using  the  radiation  centers  extracted  at  only  a  limited  number  of 
angles.  In  the  second  example,  we  carry  out  a  study  to  reduce 
platform  effects.  It  is  shown  that  the  extracted  radiation  center 
information  from  the  SBR  simulation  engine  can  be  exploited 
to  identify  those  parts  of  the  target  that  give  rise  to  significant 
^secondary  radiation.  Consequently,  by  placing  absorbers  over  a 
limited  portions  of  the  platform,  it  is  possible  to  suppress  a  ma¬ 
jority  of  the  platform  scattering. 


001 8-926X/00S 1 0.00  ©  2000  IEEE 


OZDEMIR  el  ai :  ANTENNA  RADIATION  PATTERNS  ON  A  COMPLEX  PLATFORM 


993 


Fig.  1 .  CAD  model  of  an  airplane  with  a  small  vertical  dipole  antenna  placed 
above  the  cockpit. 

It  should  be  noted  that  the  numerical  results  used  throughout 
this  paper  will  be  based  on  the  SBR  code  Apatch,  which  does 
not  include  higher  order  scattering  phenomena  such  as  higher 
order  diffractions  and  surface/creeping  waves.  Furthermore,  the 
example  configurations  (e.g.,  dipole  above  aircraft)  considered 
in  this  paper  may  not  be  construed  as  practical.  However,  we 
believe  these  limitations  do  not  significantly  detract  from  the 
main  point  of  this  paper,  which  is  the  concept  of  radiation  center 
representation  for  modeling  platform  effects.  ;  . 

II.  Radiation  Center  Extraction 

A.  ASAR  Image  Formation  Algorithm  and  Radiation  Center 
Model  ;  v 

Before  presenting  the  radiation  center  model  for  antenna  ra¬ 
diation  data,  we  shall  first  review  the  ASAR  image  formation 
algorithm  presented  earlier  in  [1]  and  provide  some  motivations 
for  the  radiation  center  model.  To  form  an  ASAR  image  from 
antenna  radiation  data,  we  first  assume  the  antenna  is  located  at 
the  origin  as  shown  in  Fig.  3.  The  scattered  electric  field  around 
-x  direction  from  a  scattering  point  P(x0,y0,  za)  on  the  plat¬ 
form  can  be  approximated  by 

Es(u>,  4>,0)  =  A-  e-Jfc(r»+I»)  ■  e~jk°*y°  •  e-ik°6z°  (1) 

where  A  is  the  strength  of  the  scattered  signal,  ra  = 
\Ao  +  Vo  +  zo  is  the  path  traveled  by  the  radiated  signal  from 
the  antenna  to  point  P,  k  is  the  free-space  wave  number,  and 
k0  is  the  wave  number  at  the  center  frequency.  A  small-angle 
small-bandwidth  approximation  is  used  to  arrive  at  the  above 
expression.  By  setting  u  =  r  +  x  and  utilizing  the  linear 
dependence  of  the  phase  terms  in  (1)  on  frequency  and  angles, 
we  can  take  the  3-D  inverse  Fourier  transform  of  the  scattered 
electric  field  with  respect  to  k,  k0<j>,  and  k09  to  generate  a  3-D 
ASAR  image  as  follows: 

ASAR(u,  y,  z)  =  IFT3{i?*(u>,  <t>,  0)} 

=  IFT3{A •  e~ikUo  ■  e~iko4,y°  ■  e~iK6z°} 

=  A  •  6(u  -  u0)  ■  6{y  —  y0)  •  S(z  —  z„).  (2) 


To  convert  the  image  into  the  original  xyz-spa.cc  from  the 
uyz- space,  we  can  use  the  transformation  formula 


Therefore,  by  Fourier  processing  of  the  multifrequency 
multi-aspect  radiation  data,  the  point  P  on  the  platform  can  be 
mapped  to  a  peak  in  the  image  at  its  correct  spatial  location 
(x0  ,y0,z0).  To  summarize,  the  ASAR  imaging  process  consists 
of  3  steps:  1)  collect  multifrequency  multi-aspect  radiation  field 
data  Es(lj,  <f>,  9)\  2)  take  the  3-D  inverse  Fourier  transform  to 
form  the  image  in  the  (tt,  y,  z)  domain;  and  3)  use  the  u-to-x 
transformation  in  (3)  to  generate  the  final  image  in  the  desired 
(x,  y,  z)  domain. 

Although  the  ASAR  formation  methodology  outlined  above 
can  be  applied  to  general  radiation  data  from  either  measure¬ 
ment  or  simulation,  we  have  also  developed  an  algorithm  espe¬ 
cially  tailored  for  the  SBR  technique  called  the  image-domain 
formulation.  This  derivation  will  serve  as  the  basis  of  our  radia¬ 
tion  center  model.  In  applying  the  SBR  technique  to  the  antenna 
radiation  problem,  rays  are  first  shot  from  the  phase  center  of 
the  antenna  and  traced  according  to  geometrical  optics.  At  the 
exit  point  of  each  ray  before  they  leave  the  platform  altogether, 
a  ray-tube  integration  is  carried  out  to  find  the  contribution  of 
each  ray  to  the  total  radiated  field  at  various  frequencies  and  ob¬ 
servation  angles.  Therefore,  the  radiated  field  can  be  expressed 
as  a  weighted  sum  of  all  the  rays  that  have  been  shot  from  the 
antenna  .«  Ig  '  •' 

E3(w,<j>,9 )  =  53  at- e~jkUi  ■  e-jk°4’Vi  ■  e~jk°6zi  (4) 

i  rays 

where  oti  is  proportional  to  the  field  strength  at  the  exit  ray-tube 
and  is  only  weakly  dependent  on  aspect  and  frequency, 
(tfi,  yuZi)  is  the  location  of  the  last  hit  point,  and  =  ft  +  x,, 
where  ft  is  the  total  path  traveled  by  ith  ray  from  the  antenna 
to  the  last  hit  point.  Since  the  ASAR  image  is  generated  from 
multifrequency  multi-aspect  data  via  Fourier  inversion,  we  can 
interchange  the  order  of  the  inverse  Fourier  transform  and  the 
ray  summation  and  carry  out  the  inverse  Fourier  transform  for 
each  ray  in  closed  form.  The  resulting  ASAR  image  is  given  by 

ASAR(x,  y,  z)  =  53  ft  •  h(x  -  xhy  -  yi,z  -  zt)  (5) 

i  rays  • 

where  -  ...w 

h(x,  y,  z)  =  e“jfco(r+x)  sine  (Afcr) 

*  sinc(A:0  A(f>y)  •  sinc(fc0  AOz),  (6) 

The  detailed  derivation  of  this  image-domain  formula  can  be 
found  in  [1]  and  will  not  be  repeated  here.  In  the  expression, 
(xi ,  yi ,  Zi)  is  the  location  of  the  hit  point  on  the  target  for  the  ith 
ray.  (For  a  multiple-bounce  ray,  it  can  be  shown  that  its  down 
range  X{  is  delayed  by  the  additional  travel  distance  from  the 
antenna  to  the  last  hit  point,  while  its  cross  range  positions 
and  Zi  correspond  to  those  of  the  last  time  point.)  ft  is  the  ray 
amplitude  AJfc,  A </>,  A 9  are  the  half-bandwidths  in  the  k-f  <f>-  and 
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Slhe'lSoLdiat^nlfni^fT^^H^01^  85Tf?  **  simulated  frequency-aspect  radiation  data  about  the  nose  using  the  code  Apatch.  (b)  The  locations 

of  the  150  radiation  centers  extracted  from  the  3-D  ASAR  image  using  CLEAN,  (c)  Reconstructed  3-D  ASAR  image  using  the  extracted  radiation  centers. 


0-domains,  respectively.  A  key  observation  from  (5)  and  (6)  is 
that  we  can  form  the  ASAR  image  using  SBR  on  a  ray-by-ray 
basis  directly  in  the  image  domain.  Each  ray  contributes  to  a 
spot  centered  about  {xuyi,Zi)  in  the  ASAR  image,  with  the 
footprint  of  the  spot  governed  by  the  3-D  ray  spread  function 
h ,  which  plays  a  role  similar  to  the  point-spread  function  used 
in  the  radar  community.  Therefore,  we  have  shown  based  on 
SBR  that  the  ASAR  image  is  composed  of  a  collection  of  point 
radiators. 

Since  tens  of  thousands  of  rays  are  traced  in  the  SBR  process, 
this  may  imply  that  tens  of  thousands  of  point  radiators  are 
needed  to  adequately  represent  an  ASAR  image.  However,  as 
we  have  observed  from  Fig.  3(a),  the  ASAR  image  of  a  complex 
platform  is  actually  quite  sparse.  This  is  because  rays  interfere 
with  one  another  to  give  rise  to  strong  coherent  scattering  over 
only  a  small  localized  portion  of  the  target.  Therefore,  an  ASAR 
image  can  be  accurately  represented  by  only  a  limited  number 


of  radiating  centers.  We  adopt  a  point  radiator  model  which  has 
exactly  the  same  form  as  (5)  to  parameterize  the  ASAR  data 

ASAR(a:,  y,  z) «.  '  .4n;  • 

n radiation  -  -  x'CM 

. ,  centers  r.  . 

;  y  ~  2/n»  z  “  2n)  (7) 

where  h(x,y,z)  is  given  in  (6).  The  remaining  task  is  to  de¬ 
termine  (:rn,  yn ,  2rn),  the  locations  of  the  radiation  centers,  and 
An ,  their  corresponding  strengths. 

B.  Extraction  and  Reconstruction  Algorithms 

r  To  carry  out  the  parameterization  in  (7),  we  apply  the  image 
processing  algorithm  CLEAN.  The  CLEAN  algorithm  is  a  pop¬ 
ular  deconvolution  technique  in  radio  astronomy  [10]  and  has 
been  successfully  utilized  for  scattering  center  extraction  [6]. 
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rig.  3.  Path  of  the  radiated  signal  due  to  a  scattering  point  P  on  the  platform. 

t  is  a  robust  iterative  procedure  that  successively  picks  out  the 
lighest  point  in  the  image,  assumes  it  is  a  radiating  center  with 
he  corresponding  strength,  and  removes  its  point  spread  re- 
jponse  from  the  image.  At  the  nth  iteration,  if  An  is  the  strength 
)f  the  highest  point  in  the  image  with  locations  (xn,  yn,  zn),  the 
5-D  residual  image  is  found  by 

}-D  Residual  Image]  n+i  =  [3-D  Residual  Image]  n 

-  [An  ■  h{x  -xn,y-  yn, z  -  zn)\. 

(8) 

Hie  extraction  process  is  iterated  until  the  maximum  in  the 
residual  image  reaches  a  user-defined  threshold.  Typically,  the 
energy  in  the  residual  image  decreases  quickly  during  the  initial 
stages  of  the  iteration  and  tapers  off  after  reaching  a  noise  floor. 

As  an  example,  the  3-D  ASAR  image  in  Fig.  2(a)  is  param¬ 
eterized  using  the  3-D  CLEAN  algorithm.  The  ASAR  image 
:s  formed  using  the  multiffequency  multi-aspect  radiation  data 
computed  using  Apatch  for  the  dipole-airplane  configuration 
shown  in  Fig.  1.  The  center  frequency  is  10  GHz  with  0.56 
GHz  of  bandwidth  sampled  over  64  frequencies.  The  observa¬ 
tion  angles  are  centered  about  the  nose  of  the  airplane  and  range 
between  —1.67°  and  1.67°  in  azimuth  sampled  over  32  points 
and  between  -0.86°  and  0.86°  in  elevation  over  eight  points, 
[n  the  ASAR  image  formation,  we  use  only  the  scattered  field 
from  the  airplane  and  not  the  direct  radiation  from  the  antenna, 
rherefore,  only  the  platform  effect  is  imaged.  Fig.  2(b)  shows 
the  locations  of  the  extracted  radiation  centers  from  the  3-D 
ASAR  image  of  Fig.  2(a)  using  the  3-D  CLEAN  algorithm.  A 
total  of  150  radiation  centers  are  extracted  from  the  image  and 
their  locations  are  plotted  in  Fig.  2(b)  as  small  circles.  We  ob¬ 
serve  that  most  of  them  are  concentrated  around  the  z  =  -2 
m  cut  which  corresponds  to  the  top  surface  of  the  airplane.  The 
remainders  are  located  on  the  two  tails.  Fig.  2(c)  shows  the  re¬ 
constructed  ASAR  image  from  the  150  extracted  radiating  cen¬ 
ters.  The  agreement  between  the  original  image  and  the  recon¬ 
structed  one  is  quite  good  over  a  dynamic  range  of  50  dB.  Once 
we  have  extracted  the  radiation  center  model  from  the  ASAR 
image,  we  can  reconstruct  the  frequency  and  aspect  radiation 
patterns  of  the  antenna-platform  structure.  The  reconstruction 


formula  is  essentially  the  Fourier  transform  of  (7)  and  is  given 
by 

E\u,  4>16)  =  y£lAn-  e"J'fe(rn+x")  •  e~ik°*y'  ■  e~iK6z- . 

n 

(9) 

By  applying  the  above  formula,  we  reconstruct  the  platform  ra¬ 
diation  patterns  in  frequency  and  aspect  using  the  150  radia¬ 
tion  centers.  Fig.  4(a)  shows  the  comparison  between  the  re¬ 
constructed  frequency  domain  data  and  the  original  Apatch  data 
calculated  over  the  0.56-GHz  bandwidth  about  10  GHz.  The  ob¬ 
servation  angle  is  along  the  nose.  Fig.  4(b)  shows  the  compar¬ 
ison  between  the  reconstructed  azimuth  data  calculated  over  the 
3.33°  window  about  the  nose.  The  frequency  is  10  GHz  and  the 
elevation  angle  is  0°.  As  shown  from  the  graphs,  good  agree¬ 
ment  is  achieved  between  the  reconstructed  data  and  the  orig¬ 
inal  data  over  the  frequency  and  aspect  from  which  the  radiation 
center  model  was  generated. 

m.  Radiation  Center  Applications 

In  the  previous  section,  it  is  shown  that  the  radiation  center 
model  for  an  arbitrary  antenna-platform  configuration  can  be 
extracted  using  the  ASAR-CLEAN  methodology.  In  this  sec¬ 
tion,  we  present  two  examples  of  the  utilities  of  such  radiation 
center  representation. 

A.  Radiation  Data  Compression 

The  antenna  radiation  pattern  on  an  electrically  large  plat¬ 
form  is,  in  general,  a  very  rapidly  varying  function  of  frequency 
and  aspect.  In  this  example,  we  shall  demonstrate  how  the  ra¬ 
diation  center  concept  developed  in  the  last  section  can  be  uti¬ 
lized  to  create  a  sparse  model  of  the  antenna-platform  radiation 
over  frequencies  and  aspect  angles,  thus  achieving  data  com¬ 
pression.  Note  that  for  a  complex  platform,  the  radiation  center 
model  extracted  at  a  particular  observation  angle  is  not  expected 
to  be  valid  over  large  angular  extent.  This  is  due  to  shadowing 
and  other  complex  multiple  scattering  phenomena.  Therefore,  to 
fully  characterize  the  antenna-platform  radiation  at  all  aspects, 
we  need  to  extract  3-D  radiation  center  models  at  various  angles 
on  a  uniform  grid  in  both  elevation  and  azimuth.  Once  the  ra¬ 
diation  center  models  are  extracted  at  all  the  angles  on  the  grid, 
we  can  obtain  the  radiated  field  at  any  arbitrary  angle  by  table 
lookup  and  reconstruction.  The  key  question  in  this  construct 
is  the  granularity  of  the  grid,  which,  in  practice,  will  depend 
on  the  platform  complexity  as  well  as  the  needed  accuracy  in 
the  reconstruction.  In  our  example,  we  restrict  our  attention  to 
the  azimuth  direction,  as  shown  in  Fig.  5.  The  antenna-platform 
configuration  is  the  same  that  shown  in  Fig.  1.  First,  the  3-D 
ASAR  images  for  different  observation  angles  in  the  azimuth 
direction  are  generated  using  Apatch.  The  angular  granularity 
is  chosen  to  be  5°  so  that  a  total  of  72  images  are  generated  in 
covering  the  full  360°  azimuth  in  the  zero  elevation  plane.  No¬ 
de^  that  only  a  one-time  ray  shoot  is  needed  from  the  antenna  in 
constructing  the  72  ASAR  images  because  of  the  bistatic  nature 
of  the  ASAR  scenario.  This  is  contrary  to  the  monostatic  radar 
signature  studies  we  have  carried  out  previously  [12]:  We  use  a 
center  frequency  of- 10  GHz  and  a  bandwidth  of  0.49  GHz.  To 


f 


996 


IEEE  TRANSACTIONS  ON  ANTENNAS  AND  PROPAGATION,  VOL.  48,  NO.  6,  JUNE  2000 


0 

*o 

p-IO 

o 

*-20 


a -30 

-40 

9.75  9.8  9.85  9.9  9.95  10  10.05  10.1  10.15  10.2  10.25 

- >  Frequency  (GHz) 

(a) 

0 

s 

•a 

^  “’10 
© 

'f 

*  -20 

1 

a  -30 

-40 

-1-5  -1  -0.5  0  0.5  1  1.5 

— >  Azimuth  (Deg.) 

(b) 

Fig.  4.  Comparison  of  the  original  and  the  reconstructed  radiation  patterns,  (a)  Frequency  sweep,  (b)  Azimuth  sweep. 


Fig.  5.  Global  radiation  center  representation  of  the  radiation  pattern  along  the  azimuth  cut  obtained  via  ASAR  image  formation  and  radiation  center  extraction 
at  various  angles.  ... 


generate  each  ASAR  image,  the  azimuth  bandwidth  is  chosen 
to  be  2.22°  and  the  elevation  bandwidth  is  taken  as  2.46°. 

Next,  we  extract  100  radiation  centers  from  each  of  the  72 
3-D  ASAR  images  using  CLEAN.  Shown  in  Fig.  6  are  four 
reconstructed  ASAR  images  at  0°,  30°,  60°,  and  90°  with  re¬ 
spect  to  the  nose  of  the  airplane.  They  are  2-D  projected  images 


formed  by  summing  up  all  the  2-slices  of  the  original  3-D  im- 
p  ages.  We  observe  that  there  are  some  persistent  hot  spots  on  the 
nose  and  the  right  wing  of  the  airplane,  which  appear  to  be  vis¬ 
ible  in  all  four  images,  while  other  features  have  smaller  visi¬ 
bility  extent  To  test  the  validity  of  the  model*  the  platform  ra-: 
diation  pattern  in  the  0°-180°  azimuth  range  is  reconstructed; 


pver  each  5°  sector  by  using  the  extracted  radiation  centers  and 
9).  It  is  shown  as  the  solid  line  in  Fig.  7(a).  Also  plotted  in 
flashed  line  in  the  figure  is  the  original  azimuth  radiation  pat¬ 
tern  computed  using  Apatch  at  a  granularity  of  1° .  Note  that  the 
reconstructed  data  can  be  computed  to  any  desired  fine  angular 
granularity  very  rapidly  from  the  radiation  center  model.  The 
Apatch  data,  on  the  other  hand,  are  computed  very  coarsely  as  it 
[s  very  time-consuming  to  generate.  The  two  curves  show  some 
Agreement  in  overall  trend,  although  the  exact  Apatch  data  are 
too  grossly  undersampled.  In  Fig.  7(b),  the  comparison  between 
the  two  curves  computed  on  the  same  fine  granularity  of  0.07° 
pver  the  30° -35°  sector  is  shown.  The  reconstructed  data  from 
30° -32.5°  are  reconstructed  from  one  set  of  radiation  centers 
while  those  from  32.5°-35°  are  reconstructed  from  another  set. 
kote  that  the  data  between31. 11°  and  33.89°  correspond  to  “ex¬ 
trapolated”  data  based  on  the  nearest-neighbor  radiation  center 
model.  The  overall  agreement  between  the  two  sets  of  data  are 
fairly  good,  even  in  the  extrapolated  region.  We  have  also  at¬ 
tempted  to  extend  the  extrapolation  to  a  10°  window,  and  have 
noticed  significant  degradation  in  the  outer  regions  of  extrapo- 

Ition.  Therefore,  the  5°  granularity  seems  to  be  a  good  choice 
r  this  example.  Reconstruction  fidelity  can  be  improved  by 
ing  more  radiation  centers  in  the  CLEAN  procedure  and  by  re- 
icing  the  angular  granularity.  This,  of  course,  is  at  the  expense 
model  sparsity.  Finally,  the  same  philosophy  can  also  be  ap- 
ied  to  the  elevation  direction  to  fully  exploit  the  scheme  for 
tta  compression.  The  data  compression  ratio  achievable  using 
e  3-D  radiation  center  model  can  be  roughly  estimated  as  fol- 
ws.  We  assume  the  original  radiation  data  cover  the  entire  vis¬ 


ible  sphere  and  64  frequency  samples.  The  original  complex  ra¬ 
diation  data  set,  sampled  at  a  granularity  of  0. 1 0 ,  requires  a  total 
storage  space  of  (360/0.1)  *  (180/0.1)  *  (64)  *  (8  bytes)  =  3.3 
Gb.  The  radiation  center  model  contains  100  radiation  centers 
per  angle,  with  each  radiation  center  requiring  the  storage  of  five 
real  numbers  (three  for  the  location  and  two  for  the  complex 
amplitude).  By  using  a  5°  granularity  to  cover  the  entire  visible 
sphere,  the  storage  space  required  for  the  radiation  center  set 
amounts  to  (360/5)  *  (180/5)  *  (100)  *  (5)  *  (4  bytes)  =  5.2 
Mb.  The  data  compression  ratio  is  about  640 : 1 .  To  summarize, 
we  have  shown  that  a  global  sparse  radiation  center  model  can 
be  constructed.  Once  available,  it  can  be  used  to  reconstruct  ra¬ 
diation  data  to  a  very  fine  granularity  with  good  fidelity, 

B.  Platform  Effect  Reduction  /t ... 

In  the  second  example,  we  shall  utilize  the  radiation  center: 
representation  to  carry  out  a  platform  effect  reduction  study. 
Since  the  radiation  centers  pinpoint  the  locations  of  the  dom¬ 
inant  secondary  scattering  on  the  platform,  one  way  to  mitigate 
such  undesirable  effects  is  to  place  absorbers  at  those  locations. 
This  is  analogous  to  carrying  out  signature  reduction  work  using 
ISAR  images  as  a  guide.  We  use  the  same  example  as  that  in 
Section  IH-A  and  first  generate  the  radiation  center  representa¬ 
tions  of  the  scattered  field  over  360°  azimuth  at  5°  angular  gran¬ 
ularity.  Next,  those  radiation  centers  that  contribute  more  than 
-25  dB  to  the  radiated  field  are  pinpointed.  To  eliminate  these 
radiation  centers,  we  identify  the  facets  on  the  target  that  corre¬ 
spond  to  the  locations  of  the  radiation  centers  and  convert  them 
from  perfect  conductors  into, perfect  absorbers.  By  using  this: 

Reproduced  From 
Best  Available  Copy 


998 


IEEE  TRANSACTIONS  ON  ANTENNAS  AND  PROPAGATION,  VOL.  48,  NO.  6,  JUNE  2000 


“‘‘“fT  P!frnlf.'h)®  *e  original  data  and  the  reconstructed  data  based  on  radiation  centers,  (a)  Coarse  comparison  between 
L  42  gSuS  ^  31  reconstructed  data  are  sampled  at  0.07°.  (b)  Detailed  comparison  between  30‘  and  35°  a,  tTeT^e 


logic,  we  end  up  with  387  absorbing  facets  on  the  platform  out 
of  approximately  8000  facets  that  describe  the  airplane.  The  new 
CAD  file  is  shown  in  Fig.  8(a),  where  the  absorbing  facets  are 
marked  in  a  darker  shade.  Fig.  9(a)  shows  the  azimuth  radiation 
pattern  from  the  platform  before  and  after  the  absorber  coating. 
The  patterns  are  computed  by  Apatch.  We  observe  that  the  plat¬ 
form  scattering  is  significantly  reduced  by  using  the  absorbers. 
However,  the  resulting  platform  radiation  level  is  not  suppressed 
below  the  -25-dB  design  level  as  expected.  In  particular  the 
platform  scattering  near  the  20°  and  170°  directions  remains 
quite  high.  This  can  be  explained  by  the  fact  that  the  ASAR 
imaging  algorithm  is  based  on  the  single-bounce  assumption. 
Multiple-bounce  mechanisms  are  not  correctly  mapped  to  the 
platform  in  the  ASAR  image,  but  are  rather  delayed  in  the  down- 
range  direction  [1]. 

To  circumvent  this  problem  and  further  reduce  the  platform 
scattering  from  multiple  bounce  mechanisms,  we  devise  a 
scheme  to  tie  each  radiation  center  back  to  the  ray  mechanism 
that  gave  rise  to  it.  The  basic  idea  is  to  save  the  hit  point 
information  of  the  ray  which  contributes  the  largest  amount  of 
energy  to  a  particular  ASAR  image  pixel  during  the  SBR-based 
image  formation  process  in  (5).  Consequently,  during  the 
CLEAN  procedure,  each  extracted  radiation  center  has  an  as¬ 
sociated  “hit  point  list”  to  allow  tie-back  of  the  radiation  center 
to  the  specific  hit  point  locations  on  the  platform.  Such  tie-back, 
information  is  only  approximate,  but  has  been  found  to  be  very 
useful  in  radar  signature  and  target  identification  applications 
[12],  It  is  important  to  point  out  that  this  information  can  be 
extracted  only  if  an  SBR-based  approach  is  used  to  generate  • 
the  ASAR  image.  Using  this  new  tie-back  information  for  each 
radiation  center,  we  can  now  place  absorbers  on  those  facets 
associated  with  the  hit  points  on  the  target  In  fact,  we  only  need 
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Fig.  8.  CAD  models  with  the  absorbing  facets  shown  in  a  darker  shade:  (a) ' 
Absorber  placement  using  the  first  scheme  in  which  the  facet  closest  to  a 
radiation  center  is  changed  into  an  absorber,  (b)  Absorber  placement  using  the 
second  scheme  in  which  the  hit  point  list  for  a  radiation  center  is  found  via 
SBR  and  the  facet  from  the  first  hit  is  changed  into  an  absorber. ~  "X ..  : 
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Fig.  9.  Azimuth  radiation  patterns  before  and  after  putting  absorbers  on  the  platform,  (a)  Absorber  placement  according  to  the  first  scheme,  (b)  Absorber  placement 
according  to  the  second  scheme. 
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Fig.  10.  3-D  AS AR  images  observed  along  the  nose  sector,  (a)  No  absorbers,  (b)  With  absorbers. 

to  put  an  absorber  on  any  one  of  the  hit  points  to  eliminate  the  the  first  scheme.  The  scattered  field  from  the  platform  comes 

contribution.  We  have  chosen  to  do  so  for  all  the  first-bounce  fairly  close  to  meeting  the  — 25-dB  design  level  at  most  of  the 

facets.  A  total  of  403  absorbing  facets  on  the  platform  are  angles.  In  particular  the  platform  scattering  at  20°  and  170° 

selected  in  this  manner  out  of  the  8000  facets.  Fig.  8(b)  shows  directions  is  well  suppressed.  Fig.  10(a)  and  (b)  shows  the 

the  absorber-coated  airplane  model.  Note  the  difference  in  the  3-D  ASAR  images  of  the  platform  scattering  along  the  nose 

absorbing  facet  locations  in  Fig.  8(a)  and  (b).  Fig.  9(b)  shows  sector  before  and  after  the  absorber  treatment,  respectively, 

the  platform  scattering  pattern  resulting  from  this  approach.  We  observe  that  the  platform  image  nearly  vanishes  after  the 

The  scattering  pattern  of  the  uncoated  platform  is  again  plotted  absorber  coating,  with  the  remaining  stray  radiation  coming 

in  the  figure  for  reference.  We  can  see  that  the  radiation  from  from  the  nose  portion  of  the  airplane.  Note  that  in  this  example 

the  platform  is  now  much  better  suppressed  when  compared  to  the  absorber  coating  is  placed  on  only  one  facet  for  each 
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radiation  center  selected,  no  matter  how  small  the  facet.  Since 
a  radiation  center  is  formed  by  the  coherent  radiation  from  the 
induced  current  over  some  finite  region  of  the  platform,  some  * 
minimum  area  may  be  required  to  more  effectively  eliminate 
the  radiation  from  each  scattering  center.  We  believe  this  is 
the  cause  for  the  remaining  stray  radiation,  especially  near  the 
nose  and  the  back  of  the  airplane,  where  the  facets  in  the  CAD 
model  are  very  small. 

In  this  example,  we  have  carried  out  an  exercise  to  suppress 
platform  radiation  by  selectively  coating  the  platform  from 
the  radiation  center  information.  The  first  scheme  is  simply 
based  on  placing  absorbers  at  the  radiation  center  location.  It 
works  well  for  suppressing  the  single-bounce  mechanisms.  The 
second  improved  scheme  is  based  on  the  hit  point  list  generated 
during  the  radiation  center  extraction  procedure.  It  improves 
performance  by  more  effectively  tracking  the  locations  of 
the  multibounce  mechanisms.  However,  this  scheme  is  more 
restrictive  since  it  requires  the  SBR  simulation  algorithm  be 
used  for  ASAR  image  formation. 

IV.  Conclusion 

In  this  paper,  we  presented  a  sparse  representation  of  antenna 
radiation  pattern  on  a  complex  platform.  This  representation  is 
based  on  a  point  radiator  model  that  describes  the  radiation  pat¬ 
tern  by  a  collection  of  radiation  centers  on  the  platform.  The 
methodology  for  obtaining  the  radiation  center  model  is  to  first 
generate  the  3-D  antenna  synthetic  aperture  radar  imagery  of 
the  platform  and  then  to  parameterize  the  resulting  image  by  a 
collection  of  point  radiators  via  the  CLEAN  algorithm.  It  was 
shown  that  once  such  a  representation  is  obtained,  we  can  re¬ 
construct  and  extrapolate  antenna  radiation  patterns  over  fre¬ 
quencies  and  aspects  with  good  fidelity,  thus  achieving  high 
data  compression  ratio.  Furthermore,  it  was  shown  that  when 
coupled  with  the  SBR  simulation  engine,  the  resulting  radia¬ 
tion  center  information  can  be  used  to  pinpoint  cause-and-ef- 
fect  in  platform  scattering  and  provide  important  guidelines  for 
reducing  platform  effects.  Finally,  we  should  point  out  that  the 
data  used  throughout  this  paper  were  based  on  ray-tracing  sim¬ 
ulation.  Validation  from  measurement  or  more  rigorous  numer¬ 
ical  data  was  not  attempted.  However,  we  have  since  carried 
out  the  same  radiation  center  modeling  using  the  computational 
electromagnetics  data  from  a  full- wave  simulator.  That  result  is 
being  reported  in  [13].  It  is  shown  that  the  conclusions  reached 
in  this  paper  holds  true  even  for  more  rigorously  computed  data. 
Measurement  results  of  ASAR  images  have  also  been  reported 
recently  in  [] 
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ABSTRACT:  Using  full-wave  computational  electromagnetics  data,  we 
demonstrate  that  the  radiated  field  from  an  antenna  on  a  complex 
platform  can  be  sparsely  represented  by  a  radiation  center  model.  An 
extraction  scheme  based  on  the  matching-pursuit  algorithm  is  devised  to 
extract  the  model  parameters  from  multifrequency,  multiaspect  radiation 
data.  To  accelerate  the  parameterization  procedure,  the  previously  devel¬ 
oped  antenna  synthetic  aperture  radar  imaging  algorithm  is  used  to 
provide  an  initial  estimate  of  the  radiation  center  position  during  the 
search.  It  is  shown  that  the  resulting  radiation  center  model  provides  a 
sparse,  physical  representation  of  the  antenna-platform  interaction. 

©  2000  John  Wiley  &  Sons,  Inc.  Microwave  Opt  Technol  Lett  26: 

4-7,  2000. 

Key  words:  radiation  center  model;  antenna-platform  interaction ; 
matching-pursuit  algorithm 

I.  INTRODUCTION 

In  antenna  analysis  and  design,  the  mounting  platform  can 
have  a  significant  effect  on  the  antenna  radiation  characteris¬ 
tics.  In  addition  to  the  forward  problem  of  characterizing  the 
antenna  radiation  in  the  presence  of  the  platform,  also  of 
importance  is  the  development  of  a  sparse  model  to  describe 
the  radiation  physics  in  antenna-platform  interaction.  To¬ 
ward  this  end,  we  have  previously  proposed  a  point  radiator 
model  to  represent  the  complex  radiation  patterns  resulting 
from  antenna-platform  interactions  [1,  2].  In  this  model,  we 
assume  that  the  platform  scattering  due  to  the  primary  an¬ 
tenna  radiation  can  be  approximately  described  by  a  set  of 
“radiation  centers”  on  the  platform.  This  concept  is  similar  to 
the  scattering  center  model  widely  used  in  radar  scattering, 
where  the  high-frequency  scattering  from  a  complex  target 
often  can  be  modeled  efficiently  by  a  sparse  set  of  point 
scatterers  on  the  target  [3-6].  To  extract  the  radiation  center 
model  from  a  multifrequency,  multiaspect  antenna  radiation 
pattern,  we  have  proposed  a  parameterization  scheme  based 
on  microwave  imaging.  The  procedure  entails  first  generating 
the  3-D  antenna  synthetic  aperture  radar  (ASAR)  imagery  of 
the  platform,  and  then  parameterizing  the  resulting  image  by 
a  collection  of  point  radiators  via  the  CLEAN  algorithm.  It 
was  shown  that,  once  such  a  representation  is  obtained,  we 
can  rapidly  reconstruct  antenna  radiation  patterns  over  fre¬ 
quencies  and  aspects  with  good  fidelity.  Thus,  such  a  model 
can  be  used  for  the  real-time  reconstruction  of  complex 
antenna  patterns  in  high-level  system  simulations.  Further¬ 
more,  the  resulting  radiation  center  information  can  be  used 
to  pinpoint  cause-and-effect  in  platform  scattering,  and  to 
provide  design  guidelines  for  antenna  placement  and  opti¬ 
mization.  However,  the  method  was  based  on  a  Fourier-based 
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algorithm  that  relied  on  a  number  of  small-angle,  small-band¬ 
width  approximations.  When  these  conditions  are  not  met, 
the  resulting  parameterization  is  no  longer  sparse.  In  addi¬ 
tion,  the  concept  was  only  tested  using  data  from  high- 
frequency  ray-tracing  simulation  [7]. 

In  this  paper,  we  overcome  the  above  deficiencies  by:  1) 
developing  a  more  generalized  algorithm  based  on  matching 
pursuit  [8]  to  extract  the  radiation  center  model  from  fre¬ 
quency-aspect  data,  and  2)  demonstrating  the  algorithm  using 
computational  electromagnetics  (CEM)  data  from  the  full- 
wave  numerical  solver  FISC  [9].  In  the  matching-pursuit 
algorithm,  the  radiation  center  model  is  extracted  via  an 
iterative  projection  process.  During  each  stage  of  the  itera¬ 
tion,  the  basis  function  associated  with  each  possible  radia¬ 
tion  center  position  is  projected  onto  the  radiation  data.  The 
point  that  maximizes  the  projection  coefficient  is  selected  as 
a  radiation  center,  and  its  contribution  is  subtracted  from  the 
total  radiated  field.  This  process  is  repeated  until  the  energy 
in  the  residual  signal  has  reached  a  sufficiently  small  level.  To 
accelerate  the  exhaustive  search  process  during  each  itera¬ 
tion,  the  ASAR  image  is  used  to  obtain  an  initial  estimate  of 
the  radiation  center  position  before  the  fine  search  for  its 
precise  location.  By  applying  this  algorithm  to  full-wave  CEM 
simulation  data,  we  demonstrate  that  the  secondary  radiation 
from  a  complex  platform  can  be  compactly  represented  by  a 
relatively  small  number  of  radiation  centers. 

II.  RADIATION  CENTER  MODEL 
AND  ITS  PARAMETERIZATION 

The  model  that  we  will  adopt  assumes  that  the  secondary 
radiation  from  the  platform  due  to  the  primary  radiation 
from  an  antenna  can  be  described  by  a  set  of  radiation 
centers.  For  each  radiation  center  located  at  position 
(*o>  yo>  (see  1),  the  corresponding  radiated  field  as  a 
function  of  frequency  and  aspect  can  be  written  as 

Er(ft  0j,  02)  =  Ae~^kr°e^x°co&  01  cos  sin  e2+z 0  sin  e,> 

(1) 

where  k  =  lirf/c  is  the  wavenumber,  r0  =  yjx\  +  yl  +  Zq  *s 
the  distance  from  the  antenna  to  the  radiation  center,  and 
the  angles  0,  and  02  are  defined  with  respect  to  the  local 
coordinate  system  in  which  the  x-axis  points  toward  the 
central  observation  direction  (see  Fig.  1).  Note  that  the  phase 
factor  in  (1)  accounts  for  the  path  delay  from  the  antenna  to 
the  radiation  center,  and  then  to  the  observation  direction  in 
the  far  field.  When  this  model  is  used  to  describe  the  sec¬ 
ondary  radiation  from  a  complex  platform,  the  total  field  can 
be  expressed  as  a  sum  of  the  fields  due  to  all  of  the  radiation 

Observation  Radiation 


Figure  1  Radiation  center  model  for  the  secondary  platform  radia¬ 
tion 


centers: 

££*.,</.  LE[(f,eue2) 

i 

=  Y,'4iQi(.f,0u62,xi,yi,zi).  (2) 

i 

The  function  <J>;  includes  the  exponential  terms  in  (1),  and  is 
the  basis  function  due  to  the  ith  radiation  center.  Given  the 
radiated  field  as  a  function  of  frequency  and  aspect,  our  goal 
is  to  determine  a  model  that  best  approximates  the  radiated 
field  with  as  few  radiation  centers  as  possible.  The  unknown 
parameters  are  the  position  of  the  radiation  centers  (xi9yi9  zt) 
and  their  strength  Ar 

Note  that,  due  to  the  presence  of  the  e~jkr°  term,  the 
usual  Fourier  relationship  does  not  exist  between  the  fre¬ 
quency-aspect  variables  and  the  position  of  the  radiation 
centers.  In  addition,  since  the  basis  functions  are  not  strictly 
orthogonal,  the  simultaneous  projection  of  the  radiation  data 
onto  these  bases  will  result  in  strong  interference  among  the 
bases.  To  overcome  this  problem,  we  adopt  the  matching- 
pursuit  algorithm  [8]  to  extract  the  model  parameters  itera¬ 
tively.  In  each  step  of  the  iteration,  we  project  the  radiated 
field  onto  every  possible  radiation  center  basis  in  3-D  space. 
The  location  that  gives  the  largest  projection  is  chosen  as  the 
strongest  radiation  center.  The  search  can  be  expressed  as 

A—  max  {<£r, <!>(*, y,z)>)  (3) 

(x,y,z) 

where  the  inner  product  is  defined  as 

<£'\<D>=  j  j  j  Eri>*dfd91de2.  .  (4) 

JfJd\  Je2 

Once  the  strongest  radiation  center  is  determined,  its  associ¬ 
ated  radiated  field  is  subtracted  from  the  total  radiated  field 
to  generate 

Ef+t-Ef-A,*,.  (5) 

Since  the  strongest  radiation  center  is  removed  from  the 
signal,  its  interference  on  the  weaker  radiation  centers  is 
minimized.  This  process  is  repeated  to  search  for  subsequent 
radiation  centers,  and  the  process  terminates  when  the  re¬ 
mainder  of  the  energy  in  the  radiated  field  has  reached  a 
sufficiently  small  level.  At  the  end  of  the  iteration  process,  a 
radiation  center  representation  is  obtained  to  represent  the 
original  radiated  field. 

A  problem  with  the  matching-pursuit  algorithm  is  that  it 
requires  exhaustive  search  over  3-D  space  during  each  itera¬ 
tion  of  the  extraction  procedure.  This  can  be  prohibitively 
time  consuming.  To  reduce  the  search  time,  we  apply  the 
approximate  Fourier-based  ASAR  algorithm  [1]  to  first  esti¬ 
mate  the  position  of  the  strongest  radiation  center  before 
carrying  out  an  exhaustive  fine  search  around  this  neighbor¬ 
hood.  The  ASAR  imaging  algorithm  is  derived  by  applying 
the  small-bandwidth,  small-aspect  approximations  to  the 
phase  terms  in  (1): 

kx  =  k  cos  cos  d2  ~  k 

ky  =  k  cos  01  sin  $2  *  kc  62 

kz  =  k  sin  0l  =  kc6x  (6) 


MICROWAVE  AND  OPTICAL  TECHNOLOGY  LETTERS  /  Vol.  26,  No.  1 ,  July  5  2000  5 


where  kc  is  the  wavenumber  at  the  center  frequency.  Substi¬ 
tuting  (6)  into  (1),  we  can  write  the  radiated  field  as 


£r(/,  0lf  02)  =Ae~Jkir«-x»)ejkeiyaejkd'z».  (7) 

Note  that  (7)  implies  a  Fourier  relationship  between  (r0  - 
xn,y0,  z0)  and  (/,  9X,  02).  Therefore,  by  Fourier  transforming 
the  radiated  field,  we  obtain  an  image  in  the  (r  —  x,  y,  z) 
space.  The  point  with  the  highest  value  in  the  image  is 
associated  with  the  strongest  radiation  center.  We  can  fur¬ 
ther  obtain  its  (x,  y,  z)  position  by  using  the  transformation 


*-2 


l((y2  +  z2) 
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(8) 
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Figure  3  Extracted  point  scatterers  for  horizontal  dipole  antenna 
on  the  platform 


where  u  =  r  -  x.  We  found  that,  when  the  small-bandwidth, 
small-angle  conditions  are  not  satisfied,  the  ASAR  image  still 
provided  a  fair  estimate  of  the  strongest  radiation  center 
position.  Therefore,  during  the  search  process  in  the  match¬ 
ing  pursuit,  we  use  the  ASAR  algorithm  to  carry  out  a  fast, 
coarse  search  of  the  position  of  the  strongest  radiation  cen¬ 
ter.  We  then  zoom  in  on  the  precise  location  of  the  radiation 
center  via  a  fine  search.  As  a  result,  the  total  number  of 
search  points  is  reduced,  and  the  resolution  can  be  improved. 

III.  NUMERICAL  RESULTS 

In  this  section,  we  validate  the  radiation  center  model  with 
antenna -platform  radiation  data  generated  via  CEM  simula¬ 
tion.  The  computation  is  carried  out  with  FISC  [8],  which  is  a 
method-of-moments  solver  based  on  the  multilevel  fast  multi¬ 
pole  method.  The  layout  of  the  antenna-platform  structure  is 
shown  in  Figure  2.  The  ship-like  platform  is  10  m  in  length 
and  2.5  m  in  height.  The  antenna  is  a  short  dipole,  and  is 
placed  2  m  above  the  deck.  We  compute  the  radiated  field 
with  FISC  at  a  center  frequency  of  1.0  GHz  and  a  bandwidth 
of  500  MHz.  The  center  observation  angle  is  <£el  —  40°, 
<t>A z  =  50°,  and  the  angular  range  is  about  23°  in  both  azimuth 
and  elevation  angles.  To  reduce  the  error  in  the  ASAR 
algorithm  resulting  from  applying  the  small-angle  approxima¬ 
tion  to  a  fairly  large  observation  range  (23°),  we  use  the 
actual  wavenumbers  (kx,  kyy  kz)  to  replace  the  frequency  and 
angles  (k,  0l9  02 )  in  (7).  The  radiated  field  is  computed  uni¬ 
formly  in  the  (kx,  ky,  kz)  space. 

We  first  consider  the  source  dipole  parallel  to  the  deck, 
and  extract  the  radiation  centers  using  the  procedure  de¬ 
scribed  in  the  previous  section.  The  first  20  radiation  centers 
are  shown  in  Figure  3.  The  strengths  of  the  radiation  centers 
are  represented  by  different  colors.  We  observe  that  the 


dominant  platform  scattering  comes  from  the  bow  of  the  ship 
and  the  top-hat  structure  formed  by  the  cylinder  and  deck. 
Once  the  radiation  center  model  is  obtained,  the  radiated 
field  easily  can  be  reconstructed  using  the  model.  Figure  4(a) 
is  a  plot  of  the  original  radiated  field  as  a  function  of  kx  and 
ky  which,  according  to  (6),  are  approximately  proportional  to 
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Figure  4  Comparison  of  original  and  reconstructed  radiation  pat¬ 
tern  at  kz  =  0.  (a)  Computed  by  FISC,  (b)  Reconstructed  from 
radiation  center  model 
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Figure  5  Extracted  radiation  centers  for  vertical  dipole  antenna  on 
the  platform 


the  frequency  and  relative  azimuth  angle,  respectively.  kz  is 
set  to  be  0.  The  field  reconstructed  from  the  first  20  radiation 
centers  is  shown  in  Figure  4(b).  We  observe  that  the  qualita¬ 
tive  agreement  between  Figure  4(a)  and  (b)  is  good.  We  also 
compute  the  correlation  index  between  the  two  figures  using 


Figure  6  Comparison  of  original  and  reconstructed  radiation  pat¬ 
tern  at  kx  =  0.  (a)  Computed  by  FISC,  (b)  Reconstructed  from 
radiation  center  model 


the  formula 


f  [Ef(kx,ky)E2(kx,ky)dkxdky 

R  —  ~r - — - •  (9) 

-j  j (!£,(*„  ky)\2  +  \E2{kx,ky)\2)  dkx  dky 

The  resulting  correlation  index  is  0.958,  showing  quantita¬ 
tively  that  the  actual  field  can  be  fairly  well  reconstructed  by 
the  radiation  center  model. 

Next,  we  consider  the  dipole  vertical  to  the  deck,  and 
repeat  the  same  procedure.  The  first  20  radiation  centers  are 
plotted  in  Figure  5.  Their  positions  indicate  that  a  strong 
interaction  still  comes  from  the  edges  and  corners  on  the 
platform.  But  compared  to  the  horizontal  dipole  example, 
there  are  no  radiation  centers  on  the  deck  directly  under  the 
antenna.  This  is  because  the  radiated  field  from  the  vertical 
dipole  is  small  in  that  direction.  The  original  radiated  field  is 
plotted  in  Figure  6(a)  as  a  function  of  ky  and  kz ,  which  are 
approximately  proportional  to  the  relative  azimuth  and  eleva¬ 
tion  angles,  respectively.  kx  is  set  to  be  0.  The  field  recon¬ 
structed  from  the  model  is  shown  in  Figure  6(b).  Again,  the 
agreement  between  the  two  is  good.  The  correlation  index 
between  the  two  is  0.978. 

IV.  CONCLUSIONS 

In  this  paper,  we  have  demonstrated,  using  full-wave  CEM 
simulation  data,  that  the  radiated  field  from  an  antenna  on  a 
complex  platform  can  be  sparsely  represented  by  a  radiation 
center  model.  A  matching-pursuit  algorithm  is  devised  to 
determine  the  model  parameters.  The  previously  developed 
ASAR  imaging  algorithm  is  used  to  accelerate  the  parameter¬ 
ization  procedure.  The  resulting  radiation  center  model  pro¬ 
vides  a  sparse,  physical  representation  of  the  antenna-plat¬ 
form  interaction. 
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Abstract— A  frequency-aspect  extrapolation  algorithm  is  pro¬ 
posed  to  accelerate  ISAR  image  simulation  using  fast  multipole 
solvers.  A  two-dimensional  (2-D)  multiple-arrival  model  based  on 
high-frequency  physics  is  proposed  to  parameterize  the  induced 
currents  on  the  target.  A  2-D  estimation  of  parameters  via  rotation 
invariance  technique  (ESPRIT)  algorithm  is  developed  to  estimate 
the  model  parameters  from  a  limited  number  of  computed  data 
samples  in  frequency  and  aspect.  The  model  is  then  extrapolated 
to  other  frequencies  and  aspects  to  arrive  at  broadband,  wide-angle 
radar  cross  section  (RCS)  data  for  inverse  synthetic  aperture  radar 
(ISAR)  image  construction.  This  algorithm  is  tested  using  a  canon¬ 
ical  cylinder-plate  structure  to  evaluate  its  performance.  The  ISAR 
image  of  the  benchmark  VFY-218  airplane  at  UHF  band  is  then 
predicted  using  the  fast  multipole  solver  FISC  and  the  2-D  extrap¬ 
olation  algorithm.  The  resulting  image  compares  favorably  with 
that  obtained  from  chamber  measurement  data. 

Index  Terms — Fast  multipole  solver,  frequency-aspect  extrapo¬ 
lation,  ISAR  image  simulation,  two-dimensional  ESPRIT. 


I.  Introduction 

INVERSE  synthetic  aperture  radar  (ISAR)  imaging  is  an  im¬ 
portant  tool  for  radar  signature  diagnostic  and  target  identi¬ 
fication  [1],  [2].  In  order  to  simulate  the  ISAR  image  of  a  target, 
it  is  necessary  to  solve  the  electromagnetic  scattering  problem  at 
multiple  frequencies  and  angles.  Then,  by  performing  a  two-di¬ 
mensional  (2-D)  Fourier  transform  on  the  resulting  frequency- 
aspect  radar  cross  section  (RCS)  data,  a  2-D  ISAR  image  of  the 
target  can  be  constructed.  In  this  work,  we  consider  ISAR  image 
simulation  of  complex  targets  using  a  full- wave  numerical  tech¬ 
nique.  The  focus  of  our  attention  is  the  class  of  iterative  solvers 
based  on  the  fast  multipole  method  [3],  [4],  These  solvers  have 
much  lower  computational  complexity  than  traditional  moment 
method  (MoM)  for  scattering  problems  at  a  single  frequency 
and  a  single  observation  angle.  For  multiple  frequency-aspect 
RCS  calculations,  however,  the  solver  has  to  be  executed  re¬ 
peatedly  for  each  angle  and  frequency. 

To  reduce  the  computation  burden  incurred  by  multiple  fre¬ 
quency-aspect  calculations,  the  concept  of  model-based  param¬ 
eter  estimation  can  be  applied  to  populate  the  required  data  set 
from  a  sparse  set  of  computed  data  [5]— [12].  In  particular,  we 

Manuscript  received  August  30, 1999;  revised  April  6, 2000.  This  work  was 
supported  in  part  by  the  Air  Force  MUR1  Center  for  Computational  Electro¬ 
magnetics  under  Contract  AFOSR  F49620-96- 1-0025  and  the  Office  of  Naval 
Research  under  Contract  NOOO 14-98- 1-0 178. 

The  authors  are  with  the  Department  of  Electrical  and  Computer  En- 
gineering.  University  of  Texas,  Austin,  TX  78712-1084  USA  (e-mail: 
ling@ece.utexas.edu). 

Publisher  Item  Identifier  S  0196-2892(00)05909-X. 


have  recently  developed  a  one-dimensional  (1-D)  frequency  ex- 
trapolation  algorithm  [10],  [12]  and  a  1-D  angular  extrapolation 
algorithm  [11]  to  predict  the  high-frequency  RCS  of  complex 
targets.  Our  approach  is  based  on  modeling  the  induced  current 
on  the  target  using  a  multiple-arrival  model  that  closely  resem¬ 
bles  the  ray-optical  behavior  at  high  frequencies.  The  model 
coefficients  are  determined  by  the  superresolution  estimation 
of  parameters  via  rotation  invariance  technique  (ESPRIT)  algo¬ 
rithm  [13],  [14].  In  this  paper,  we  extend  our  algorithm  to  two 
dimensions  to  carry  out  simultaneous  frequency-aspect  extrap¬ 
olation.  We  begin  by  proposing  a  2-D  multiple-arrival  model  in 
the  frequency-aspect  domain.  Next,  we  implement  a  2-D  ES¬ 
PRIT  algorithm  [15],  [16]  to  estimate  the  model  coefficients 
from  a  limited  number  of  frequency-aspect  current  data  com¬ 
puted  using  a  fast  multipole  solver.  Once  the  model  is  deter¬ 
mined,  the  currents,  and  subsequently  the  RCS  at  other  frequen¬ 
cies  and  aspects,  can  thus  be  computed. 

This  paper  is  organized  as  follows.  In  Section  n,  we  examine 
the  frequency  and  angular  dependency  of  the  current  phase  and 
propose  a  multiple-arrival  model  in  the  frequency-aspect  do¬ 
main.  In  Section  III,  a  2-D  ESPRIT  algorithm  is  developed  and 
applied  to  estimate  the  time-of-arrival  and  cross  range  param¬ 
eters  in  the  frequency-aspect  model  from  a  few  calculated  fre¬ 
quency-aspect  samples.  In  Section  IV,  the  extrapolation  algo¬ 
rithm  is  tested  bn  a  canonical  2-D  cylinder-plate  structure.  The 
result  is  compared  against  both  the  exact  result  and  that  ob¬ 
tained  by  using  the  well-known  bistatic  approximation  in  angle 
[12],  [17],  [18].  The  extrapolated  results  are  also  generated  as 
a  function  of  the  number  of  calculated  samples  to  investigate 
the  convergence  of  the  algorithm.  Finally,  we  apply  the  algo¬ 
rithm  to  extrapolate  the  RCS  and  generate  a  2-D  ISAR  image 
of  the  VFY-218  airplane  at  UHF  band.  The  extrapolated  result 
is  compared  to  the  image  obtained  from  chamber  measurement 
data  [19].  Conclusions  are  given  in  Section  V. 

II.  Two-Dimensional  Model  for  Frequency-Aspect 
Extrapolation 

We  shall  first  formulate  a  physical  model  for  the  current  in¬ 
duced  on  the  surface  of  a  target  due  to  an  incident  wave.  The 
induced  current  at  any  point  is  in  general  the  result  of  multiply 
incident  waves  from  both  the  direct  excitation  and  multiple  scat¬ 
tering  from  other  parts  of  the  target.  Therefore,  we  postulate  that 
,the  induced  current  can  be  written  as  a  sum  of  multiple  incident 
waves,  each  with  a  different  travel  paths,  as  shown  in  Fig.  1.  If 
we  denote  the  down-range  direction  with  respect  to  the  incident 
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Mech.  3  Mech.  1  Mech.  2 


Fig.  1 .  Multiple-arrival  model  for  the  induced  current. 

wave  as  x  and  the  cross  range  direction  as  y ,  the  current  at  S  as 
a  function  of  frequency  /  can  be  written  as 

K 

J(f,S)  =  'Z/Ake-j2-¥dk,dk  =  xk  +  lk  (1) 

k= 1 


where  kx  =  27r/cos0/c  and  ky  =  2irfsm0/c.  Note  that  (3) 
can  be  further  approximated  to  completely  decouple  the  fre¬ 
quency  and  aspect  variable 

J(f,  0,S)  =  £  Ake~ji^liVkiee-j'i¥dki  (4) 
*  =  1 

if  we  use  sin  9  &  6  and  replace  the  frequency  variable  in  the 
first  exponential  by  the  center  frequency  fc.  Equation  (4)  then 
reduces  to  (3b)  in  [  1 1  ],  when  only  aspect  variation  is  considered. 
For  the  2-D  extrapolation  in  this  work,  we  choose  to  use  the 
model  in  (3),  since  it  is  more  accurate  when  the  aspect  range 
is  large.  Equation  (3)  is  in  the  form  of  a  sum-of-exponential 
model,  with  linear  phase  dependence  with  respect  to  both  kx  and 
ky.  Consequently,  the  2-D  superresolution  algorithm  ESPRIT 
[17]  can  be  applied  to  equally  spaced  kx  and  ky  data  samples. 
Next,  we  shall  describe  the  implementation  of  the  2-D  ESPRIT 
algorithm  to  estimate  the  unknown  parameters  in  the  model. 


where  K  is  the  number  of  incident  waves  arriving  at  5,  and  c 
is  the  speed  of  light  in  free  space.  In  the  above  definition  of  the 
path  length  dk>  we  let  (xk,  yk)  be  the  first  hit  point  on  the  target 
due  to  the  incident  wave  and  lk  be  the  total  intermediate  path 
length  of  the  multiple  scattering  mechanism  from  the  first  hit 
point  to  point  S.  Ak  is  the  amplitude  coefficient  for  each  mech¬ 
anism  and  is  assumed  to  be  frequency  independent.  Among  the 
three  mechanisms  illustrated  in  Fig.  1,  mechanism  1  is  the  direct 
incident  wave  from  the  source.  Therefore,  lx  =  0,  and  (xx,  yx) 
corresponds  to  point  S .  Mechanisms  2  and  3  are,  respectively,  a 
once-scattered  and  a  twice-scattered  wave  before  arriving  at  S. 
At  high  frequencies,  this  model  is  expected  to  be  quite  sparse, 
i.e.,  only  a  few  terms  are  needed  to  adequately  describe  the  scat¬ 
tering  physics  at  an  arbitrary  point  S  on  the  target  surface.  This 
1-D  model  was  used  to  achieve  frequency  extrapolation  at  a 
fixed  aspect  angle  in  our  earlier  work  [10],  [12]. 

Angular  dependency  can  also  be  incorporated  in  the  above 
model.  We  assume  that  all  the  intermediate  scattering  points  and 
the  amplitude  coefficient  for  each  mechanism  remain  fixed  as 
the  incident  angle  is  varied,  as  illustrated  in  Fig.  1.  This  assump¬ 
tion  was  found  to  be  fairly  accurate  for  ray-optical  fields  under 
small  angular  variation  [20].  We  have  also  applied  it  to  achieve 
angular  extrapolation  at  a  fixed  frequency  for  iterative  moment 
solvers  [11].  When  we  combine  the  aspect  behavior  with  the 
frequency  model  in  (1),  we  arrive  at  the  following  model  for  the 
current  as  a  function  of  both  frequency  and  incident  angle 

K 

J(/,0,S)  =  ^2Ake~jl 7£(*fc<cos^+yfc<sintf+/fci)#  (2) 

k=\ 

Equation  (2)  contains  three  unknowns  in  the  phase  function  of 
each  mechanism.  Next,  we  use  the  small-angle  approximation 
cos0  «1  to  arrive  at  an  expression  with  two  unknowns  in  the 
phase 

K 

J(f,  e,S)&Y,  Ake-i^dki  cos6+y *• sin9) 

k-1 

K 

=  Ake-j{kxdki+k*yki)  (3) 

k- 1 


III.  Two-Dimensional  ESPRIT  Algorithm 

In  [16],  a  2-D  ESPRIT  algorithm  was  developed  for  esti¬ 
mating  the  direction-of-arrival  in  2-D  antenna  array  problems. 
Here,  we  shall  show  that  with  some  minor  modifications,  it  can 
be  applied  to  estimate  the  parameters  Ak,  dk  and  yk  in  (3).  First, 
we  assume  that  the  parameters  are  to  be  estimated  from  known 
current  values  solved  at  M  x  N  equally  spaced  samples  in  the 
kx  -  ky  plane,  where  M  is  the  number  of  samples  in  kx  and 
N  is  the  number  of  samples  in  ky.  These  M  x  N  samples  are 
similar  to  the  elements  of  a  2-D  antenna  array  described  in  [16]. 
We  define 

pk  =  e-jAk'dk,  qk  =  e->Ak**k  (5) 

where  A kx  and  A ky  are  the  sampling  intervals  in  kx  and  fcy, 
respectively.  If  we  shift  the  origin  of  the  variables  kx  and  ky  to 
zero,  we  can  rewrite  (3)  into  a  form  similar  to  (9)  of  [16] 

K 

Zil  =YlSkPk~l(lk~l  +  nil’ 

Jk=l 

i  =  1,2,...,  M ,  /  =  1,2,...,  AT  (6) 

where  zu  is  the  current  at  the  zth  kx  value  and  the  Ith  ky  value. 
sk  is  the  modified  amplitude  coefficient  for  the  kth  mecha¬ 
nism  to  account  for  the  origin  shift,  nu  is  assumed  to  be  white 
Gaussian  noise,  which  in  our  case  is  used  to  model  the  numer¬ 
ical  error  in  the  current  computation.  Since  ESPRIT  postulates 
such  a  sum-of-exponential  model  with  additive  white  noise,  an 
averaging  procedure  has  to  be  performed  to  smooth  out  the 
noise  to  obtain  the  correct  estimate  of  the  covariance  matrix.  In 
[16],  the  averaging  is  performed  naturally  in  the  time  domain. 
However,  for  the  problem  at  hand,  the  time  dimension  does  not 
exist  and  has  to  be  synthesized.  This  can  be  accomplished  via 
a  sub-array  processing  technique  [21].  As  shown  in  Fig.  2,  a 
sjib-array  size  is  chosen  to  be  half  of  the  original  array.  Then 
shifting  the  sub-array  one  data  sample  at  a  time  in  either  kx  or 
ky  will  result  in  a  new  sub-array,  which  can  be  considered  as 
the  array  at  a  new  time  index.  If  we  assume  M  and  N  are  even 
numbers,  the  total  number  of  sub-arrays  that  can  be  generated 
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Fig.  2.  Sub-array  processing. 


in  this  manner  is  (M/2  +  1)  (N/2  4- 1).  The  signal  vectors  z (t), 
s(t)y  and  n(t)  are  defined  as  the  reshaped  column  vectors  of  the 
£th  sub-array  as  shown  in  (7)-{9),  at  the  bottom  of  the  page.  In 
the  above  expressions,  I  =  M/2  and  I  =  N/2  are  the  max¬ 
imum  row  and  column  indices  of  each  sub-array.  The  array  co- 
variance  matrix  can  thus  be  defined  exactly  like  formula  (15) 
in  [16].  Following  the  procedures  in  [16],  we  estimate  the  pa¬ 
rameters  pk  and  qk.  The  time-of-arrival  parameter  dk  and  the 
cross-range  parameter  yk  are  obtained  using 


dk 


Vk  = 


Ah 


■£qk- 


(10) 


Finally,  the  amplitude  coefficient  Ak  can  be  generated  by 
solving  (3)  via  a  least  squares  procedure.  Noted  that  the 
maximum  order  number  K  permitted  by  the  algorithm  is 


K  <  min 


<l,) 


For  the  scattering  problems  we  have  examined,  K  is  usually 
chosen  to  be  3  or  4.  The  minimum  number  of  frequency  and  an¬ 
gular  samples  can  then  be  selected  to  satisfy  the  above  criterion. 

After  the  coefficients  Ak,  dk,  and  yk  are  estimated,  we  as¬ 
sume  (3)  also  holds  for  other  frequencies  and  aspects.  There¬ 
fore,  the  induced  current  can  be  extrapolated  to  other  kx  and 
ky  values  of  interest.  An  ISAR  image  is  then  generated  using 
standard  Fourier  processing  of  the  kx  -  ky  data.  To  summa¬ 
rize,  the  extrapolation  procedure  is  earned  out  by  first  selecting 
a  number  of  densely  sampled  points  in  a  limited  frequency-as¬ 
pect  range  and  solving  the  scattering  problems  at  these  points. 
Usually,  the  computed  frequencies  are  chosen  close  to  the  low 
frequency  end,  and  the  aspect  angles  are  centered  about  the  cen¬ 
tral  angle  of  interest.  Second,  the  2-D  ESPRIT  algorithm  is  ap¬ 
plied  to  estimate  the  model  parameters  of  the  current  at  each 
point  on  the  target  surface.  Third,  the  induced  current  is  extrap¬ 
olated  to  a  wider  kx  and  ky  range  based  on  the  model  coeffi¬ 
cients  generated  from  2-D  ESPRIT  processing.  Finally,  the  cur- 


Fig.  3.  Geometry  of  the  cylinder-plate  target. 
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Fig.  4.  Reference  ISAR  image  generated  using  the  brute-force  MoM 
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Fig.  5.  ISAR  image  generated  from  71  calculated  frequency  points  and  the 
bistatic  approximation.  /  '  • 
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Fig.  6.  ISAR  images  generated  from  the  2-D  extrapolation  algorithm  using  (a)  8  x  8  calculated  points,  (b)  9  x  9  calculated  points,  (c)  10x10  calculated  points, 
and  (d)  1 1  x  11  calculated  points. 


rent  is  integrated  to  generate  the  scattered  far  field  as  a  function 
of  frequency  and  aspect  and  the  results  are  used  to  generate  the 
desired  ISAR  image. 

IV.  Numerical  Results 

To  validate  the  frequency-aspect  extrapolation  algorithm,  we 
first  consider  a  2-D  cylinder-plate  structure,  shown  in  Fig.  3. 
The  diameter  of  the  cylinder  is  4.2  m,  and  the  length  of  the 
plate  is  20  m.  The  origin  of  the  cylinder  and  the  center  of  the 
plate  are  separated  by  6.2  m.  The  frequency  band  of  interest  is 
from  0.3  GHz  to  0.65  GHz,  and  the  observation  angle  is  from 
25°  to  65°.  Fig.  4  shows  the  reference  ISAR  image  generated 
from  71  x  81  =  5751  computed  points  in  the  frequency-as¬ 
pect  plane.  The  image  has  a  dynamic  range  of  40  dB.  In  this 
structure,  strong  multiple  scattering  mechanisms  between  the 
cylinder  and  the  plate  dominate  the  backscattering.  We  can  see 
in  the  image  the  features  corresponding  to  the  direct  scattering 
from  the  cylinder,  labeled  (i),  the  front  edge  point  of  the  plate 
(ii),  and  the  shadow  boundary  cast  on  the  plate  by  the  cylinder 


(iii).  Additionally,  there  are  other  range-delayed  features  cor¬ 
responding  to  the  multiple  scattering  mechanisms.  For  com¬ 
parison,  we  calculate  the  current  at  71  frequency  samples  for 
one  aspect  angle  and  use  the  well-known  bistatic  approxima¬ 
tion  [17],  [18]  to  extrapolate  the  RCS  at  other  aspect  angles. 
Fig.  5  plots  the  resulting  ISAR  image  from  the  extrapolated 
frequency-aspect  data  using  the  bistatic  approximation.  In  this 
image,  only  the  direct  scattering  features  [(i)  and  (ii)J  are  cor¬ 
rectly  predicted,  while  the  higher-order  scattering  features  are 
poorly  predicted  in  either  position  or  amplitude.  This  is  because 
the  bistatic  approximation  is  based  on  physical  optics,  in  which 
the  current  is  assumed  to  be  excited  by  only  the  direct  incident 
wave. 

We  next  construct  the  ISAR  image  using  the  2-D  fre¬ 
quency-aspect  extrapolation  algorithm.  We  choose  8x8  points 
/at  the  low  frequency  end  and  about  the  central  angle  of  the  orig¬ 
inal  71  x  81  points.  After  the  time-of-arrival  and  cross-range 
parameters  are  extracted  by  using  2-D  ESPRIT,  the  induced 
currents  and  scattered  far  fields  are  extrapolated  to  the  kx  -  ky 
aperture  of  the  original  samples.  In  this  manner,  the  ISAR 
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Fig.  7.  IS  AR  image  of  the  VFY-2 1 8  at  1 30°  from  nose-on,  generated  from  the 
extrapolated  result  using  2-D  ESPRIT  with  36  FISC-computed  points. 


Fig.  8.  ISAR  image  of  the  VFY-2I8  at  130°  from  nose-on  generated  from 
chamber  measurement  data  [19]. 


image  can  be  generated  directly  by  a  2-D  FFT  of  the  kx  —  ky 
data  without  the  polar  reformatting  operation  [2].  The  resulting 
ISAR  image  is  plotted  in  Fig.  6(a).  As  we  can  see,  most  of  the 
features  are  predicted  in  the  correct  position.  However,  some 
of  the  weaker  features  are  highly  defocused.  We  repeat  this 
process  using  2-D  ESPRIT  extrapolation  for  9  x  9,  10  x  10, 
and  1 1  x  11  input  points.  The  resulting  images  are  shown  in 
Fig.  6(bHd),  respectively.  From  this  series  of  images,  we  see 
that  the  strong  scattering  features  are  quite  stable  throughout 
while  the  weaker  features  begin  to  converge  as  the  number  of 
input  points  is  increased.  We  find  that  the  correlation  indices 
between  the  extrapolated  images  and  the  reference  image  are, 
respectively,  0.72,  0.80,  0.82  and  0.88.  This  shows  a  steady 
improvement  as  the  number  of  input  points  is  increased.  For 
the  case  of  1 1  x  11  input  points,  nearly  all  of  the  features  in 
the  extrapolated  image  agree  well  with  those  in  the  reference 
image  from  the  brute-force  calculation.  This  corresponds  to  a 
7:1  extrapolation  ratio  in  each  dimension. 


Next,  we  use  the  2-D  extrapolation  algorithm  to  simulate  the 
ISAR  image  of  the  benchmark  VFY-2 18  airplane.  The  simula¬ 
tion  is  carried  out  using  the  multilevel  fast  multipole  code  FISC 
[22]  on  a  Pentium  II 450  MHz  PC.  The  total  computation  time 
would  take  about  200  days  if  a  brute-force  calculation  is  car¬ 
ried  out  over  the  required  41x41  frequency-aspect  samples  (to 
achieve  0.5  m  resolution  in  range  and  cross  range).  Based.on  the 
7: 1  criterion  found  above,  the  actual  computation  is  only  carried 
out  for  6  x  6  =  36  points  from  267  MHz  to  297  MHz  in  fre¬ 
quency  and  127°  to  132°  in  aspect.  The  total  computation  time 
is  40  h.  Note  that  the  time  savings  over  the  brute-force  calcula¬ 
tion  is  even  greater  than  the  (41/6)2  ratio.  This  is  because  these 
36  points  are  chosen  at  the  low  end  of  the  frequency  band  and 
take  far  less  time  to  compute  than  those  at  the  high  frequency 
end.  The  induced  current  and  far  fields  are  then  extrapolated  to  a 
kx  -  ky  grid  from  267  MHz  to  533  MHz  in  frequency  and  from 
1 10°  to  140°  in  aspect.  The  resulting  ISAR  image  of  the  airplane 
at  130°  is  shown  in  Fig.  7.  Fig.  8  shows  the  ISAR  image  con¬ 
structed  from  chamber  measurement  data  at  the  same  look  angle 
[19].  Comparing  these  two  images,  we  observe  that  nearly  all  of 
the  key  features  in  the  measurement  image  are  predicted  in  the 
simulated  image  from  using  FISC  and  the  2-D  extrapolation  al¬ 
gorithm.  As  expected,  the  image  in  Fig.  7  is  much  cleaner  than 
the  previously  predicted  image  shown  in  [12,  Fig.  11],  which 
was  simulated  from  the  1-D  frequency  extrapolation  procedure 
in  combination  with  the  bistatic  approximation.  The  improve¬ 
ment  is  due  to  the  fact  that  the  multiple  scattering  mechanisms 
are  correctly  incorporated  in  both  dimensions  in  our  new  model. 
There  still  exist  some  low  dynamic  range  noises  in  the  predicted 
image,  which  can  be  improved  with  more  computed  points. 

V.  Conclusion 

A  2-D  frequency-aspect  extrapolation  algorithm  has  been  de¬ 
veloped  to  accelerate  ISAR  image  simulation  using  fast  multi¬ 
pole  solvers.  A  2-D  multiple-arrival  model  consistent  with  high- 
frequency  scattering  physics  was  proposed  to  model  the  induced 
currents  on  the  target.  A  2-D  ESPRIT  algorithm  was  developed 
to  estimate  the  model  parameters  from  a  limited  number  of  com¬ 
puted  data  samples  in  frequency  and  aspect.  The  model  was 
then  extrapolated  to  other  frequencies  and  aspects  to  amve  at 
broadband,  wide-angle  RCS  data  for  ISAR  image  construction. 
This  algorithm  has  been  tested  using  a  canonical  cylinder-plate 
structure  to  evaluate  its  performance.  It  was  found  that  a  7:1  ex¬ 
trapolation  ratio  in  each  dimension  can  been  achieved.  Further¬ 
more,  the  input  data  samples  can  be  chosen  in  the  low  frequency 
range  for  greater  computational  payoff.  Lastly,  the  ISAR  image 
of  the  benchmark  VFY-218  airplane  at  UHF  band  has  been  pre¬ 
dicted  using  the  fast  multipole  solver  FISC  and  the  2-D  extrapo¬ 
lation  algorithm.  The  resulting  image  compared  favorably  with 
that  obtained  from  chamber  measurement  data.  The  computa¬ 
tion  time  savings  over  the  brute-force  calculation  was  about  two 
orders  of  magnitude. 

f 
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and  the  MUSIC  superresolution  algorithm.  The  c  arious  scenarios  studied 
include  element  spacing,  noise  level,  antenna  gain  and  delay,  and 
platform  effects.  Our  results  show  quantitatii  ely  the  effect  of  mutual 
coupling  on  the  direction-of-arrival  estimates.  The  compensation  of  the 
coupling  effect  using  the  coupling  matrix  approach  is  also  examined.  This 
approach  is  found  to  be  quite  satisfactory  in  most  cases,  except  when  the 
array  calibration  data  contain  a  high  noise  level  or  when  strong  platform 
effects  are  present.  ©  2000  John  Wiley  &  Sons,  Inc.  Microwave  Opt 
Technol  Lett  26:  331-336,  2000. 

Key  words:  array  mutual  coupling;  direction-of-arrival  estimates; 
compensation  by  coupling  matrix 

l  INTRODUCTION 

It  is  well  known  that  the  mutual  coupling  between  antenna 
elements  in  an  array  can  strongly  affect  their  radiation/re¬ 
ceiving  characteristics  [1,  2].  For  direction-finding  applica¬ 
tions,  the  direction-of-arrival  estimates  can  be  very  sensitive 
to  mutual  coupling,  and  such  effect  needs  to  be  properly 
accounted  for  [3].  An  effective  way  to  describe  and  compen¬ 
sate  for  the  coupling  effect  in  the  array  signal  processing 
community  is  through  the  use  of  a  coupling  matrix  [4].  It 
relates  the  active  element  patterns  of  the  individual  elements 
in  the  presence  of  the  array  environment  to  the  idealized, 
free-standing  element  patterns.  By  measuring  the  actual  ar¬ 
ray  response  at  a  few  known  incident  angles  during  calibra¬ 
tion,  the  coupling  matrix  can  be  estimated.  Such  an  array 
calibration  technique  has  been  proven  effective  in  many 
simulation  and  measurement  results  [5-7]  involving  simple 
antenna  structures. 

In  this  paper,  we  carry  out  a  study  to  simulate  the  effect  of 
mutual  coupling  for  circular  arrays.  Our  objectives  are 
twofold.  First,  we  set  out  to  simulate  the  degradation  effect 
due  to  mutual  coupling  in  direction  finding.  Second,  we  set 
out  to  examine  the  validity  of  the  coupling  matrix  model.  This 
study  is  motivated  by  anomalies  observed  in  measurement 
results  from  a  smart  antenna  testbed  for  wireless  communi¬ 
cations  [8].  Our  approach  is  to  use  the  full-wave  electro¬ 
magnetic  solver  NEC  [9]  to  carry  out  the  simulation.  The 
superresolution  algorithm  MUSIC  is  applied  to  perform  the 
direction-of-arrival  (DOA)  estimate.  The  various  scenarios 
studied  include  element  spacing,  noise  level,  antenna  gain 
and  delay,  and  platform  effects.  Our  results  show  quantita¬ 
tively  the  effect  of  mutual  coupling  on  DOA  estimates.  In 
addition,  it  is  shown  that  the  coupling  matrix  model  is  quite 
adequate  under  most  scenarios  examined  by  us. 
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II.  FORMULATIONS  FOR  DIRECTION  FINDING 

The  mutual  coupling  effect  in  an  antenna  array  is  commonly 
described  in  array  signal  processing  by  the  following  model 
[4]: 


where  Atruc  is  the  actual  array  response  matrix,  Athco  is  the 
ideal  array  response  matrix  in  the  absence  of  mutual  cou¬ 
pling,  and  C  is  an  angular-independent  coupling  matrix.  Each 
row  of  Atruc  is  the  relative  received  signal  strength  of  a 
particular  antenna  element  as  a  function  of  the  incident 
angle  of  an  incoming  plane  wave.  In  the  context  of  standard 
antenna  terminology,  this  is  known  as  the  active  element 
pattern  of  the  array  element.  In  the  absence  of  mutual 
coupling,  the  receiving  pattern  for  a  single  antenna  element 
is  independent  of  the  others,  and  is  only  a  function  of  the 
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incident  angle  and  the  position  of  the  element.  We  consider 
the  case  where  the  sources  as  well  as  the  array  are  located  in 
the  (x,  y)-plane.  Then  the  ideal  element  pattern  can  be 
written  as 

«theo ttn(<t>)  =  e^cos<t>+y-sin<t>\  m  =  1, ...»  A/  (2) 

where  M  is  the  total  number  of  elements,  (xmiym)  are  the 
coordinates  of  the  mth.  element,  and  </>  is  the  incident  angle. 
The  ideal  array  response  is  then 

^theo  “  [fltheo,  ](  $) # theo>  2^  *  (3) 

Since  Atheo  is  known  based  on  the  array  geometry,  if  Atrue  is 
given,  then  C  can  be  determined  by  using  the  pseudoinverse 
concept  as 

C  =  AtmeA^heo(AtheoA^heo)  .  (4) 

For  simple  antennas  whose  current  distributions  have  the 
same  shape  but  different  amplitudes,  C  can  be  found  from 
Atnie  at  only  a  few  observation  angles.  Once  C  is  known,  Atrue 
can  be  interpolated  to  very  fine  angular  granularity  based  on 
the  model.  Therefore,  the  coupling  matrix  concept  is  a  very 
efficient  description  of  the  mutual  coupling  effect.  We  will 
use  numerical  simulation  data  to  examine  the  validity  of  this 
model  for  direction  finding. 

To  extract  the  directions  of  arrival  from  the  received 
signals  at  the  array,  we  use  the  superresolution  algorithm 
MUSIC.  If  the  signals  come  from  N  unknown  directions,  the 
received  signal  at  the  mth  antenna  can  be  written  as 

N 

xmU)  =  £  am(4>n)sn(t )  +  zjt)  (5) 

rt*=  1 

where  sn(t )  is  the  amplitude  of  the  nth  signal,  am(<£„)  is  the 
active  element  pattern  of  the  mth  element  at  the  nth  inci¬ 
dent  angle,  and  zm(t)  represents  the  noise.  Expressed  in 
matrix  form,  (5)  becomes 

X  =  AS  +  Z  (6) 

where  A  is  the  array  response.  In  principle,  the  actual  array 
response  Atrue  should  be  used  in  (6)  to  properly  model  the 
received  signal.  However,  if  the  ideal  array  response  Atheo  is 
used  instead,  we,  in  effect,  ignore  the  mutual  coupling  be¬ 
tween  the  elements,  and  generate  what  we  will  call  the 
uncompensated  result.  In  the  next  section,  we  evaluate  the 
performance  difference  between  the  compensated  and  un¬ 
compensated  results  to  assess  the  significance  of  mutual 
coupling. 

If  the  noise  is  assumed  to  be  white  Gaussian  noise,  the 
correlation  matrix  of  X  can  be  separated  into  the  signal  part 
and  the  noise  part: 

R  =  XX"  =  ASShAh  +  an2I  (7) 

where  cr2  is  the  noise  power.  If  the  signals  are  uncorrelated, 
the  N  largest  eigenvalues  of  R  correspond  to  the  signal 
subspace,  and  the  rest  correspond  to  the  noise  subspace.  Let 
Vn  denote  the  eigenvectors  corresponding  to  the  noise  space. 


The  power  spectrum  defined  in  the  MUSIC  algorithm  is  then 

H<t>)  =  £^,^Wvnv>ra(<*>) '  (8) 

The  DOAs  are  determined  by  observing  the  peaks  in  the 
power  spectrum  as  a  function  of  incident  angle. 

111.  SIMULATION  RESULTS 

A  seven-element  circular  array  shown  in  Figure  1  is  used  in 
the  simulation.  The  radius  of  the  circle  is  0.5  wavelength.  The 
array  elements  are  half-wave  dipoles.  Each  dipole  is  divided 
into  21  segments  in  the  NEC  simulation,  with  the  load 
impedances  equal  to  73  H  at  the  center  segments.  First,  the 
actual  array  response  is  computed  densely  from  0  to  360°  in 
steps  of  0.2°.  This  will  be  used  to  generate  the  reference 
DOA  results.  Next,  we  use  the  actual  array  response  at  eight 
incident  angles  equally  spaced  around  a  circle  to  determine 
the  coupling  matrix  C  using  (4).  An  interpolated  version  of 
the  array  response  is  then  obtained  using  (1).  The  DOA 
results  generated  in  this  manner  will  be  termed  the  “com- 
pensated-by-coupling  matrix”  (CC)  results. 

To  illustrate  the  mutual  coupling  effect  in  direction  find¬ 
ing,  we  first  consider  one  incoming  signal  from  the  angle  0, 
45,  or  90°.  The  signal-to-noise  ratio  (SNR)  is  set  to  be  30  dB. 
The  normalized  MUSIC  power  spectra  are  plotted  in  Figure 
2(a).  The  dashed  curves  are  the  uncompensated  results  com¬ 
puted  using  the  ideal  array  response  Atheo.  The  solid  curves 
are  the  reference  results  generated  using  the  actual  array 
response  Atruc.  The  mutual  coupling  effect  can  be  easily 
observed  from  the  uncompensated  curves  where  the  peaks  in 
the  power  spectra  are  significantly  broadened.  It  is  interesting 
to  note  that,  due  to  the  symmetry  of  the  circular  array 
structure,  the  direction  of  arrival  is  still  estimated  correctly 
by  the  uncompensated  curve  for  a  single  incoming  signal. 
However,  as  will  be  shown  shortly,  this  will  not  be  so  when 
there  is  more  than  one  incoming  signal.  We  also  test  the 
validity  of  the  coupling  matrix  model  by  generating  the  power 
spectra  using  the  interpolated  Atrue.  The  compensated-by- 
coupling  matrix  (CC)  results  nearly  overlay  the  reference 
results,  and  are  not  shown.  Only  eight  samples  of  the  actual 
array  response  are  used  to  generate  the  7x7  coupling 
matrix.  When  compared  to  the  densely  computed  actual  array 
response  of  size  1800  X  7,  the  coupling  matrix  model  in  (1)  is 
indeed  a  very  efficient  representation  of  the  mutual  coupling 
effect.  Figure  2(b)  illustrates  the  power  spectra  when  two 
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(a) 


(b) 

Figure  2  MUSIC  power  spectrum  before  and  after  compensation  of  the  mutual  coupling  effect.  — Uncompensated, . compensated 

by  coupling  matrix  (CC),  — reference,  (a)  Power  spectrum  with  one  incoming  signal,  (b)  Power  spectrum  with  two  incoming  signals 


uncorrelated  signals  exist  at  30  and  50°.  The  dashed  curve  is 
the  uncompensated  result,  and  the  dash-dot  curve  is  the  CC 
result.  When  there  is  more  than  one  signal,  the  coupling 
effect  becomes  more  severe.  The  peaks  of  the  uncompen¬ 
sated  result  no  longer  indicate  the  correct  directions  of 
arrival.  However,  the  CC  result  still  has  sharp  peaks  at  the 
correct  DOA. 

Next,  we  examine  the  mutual  coupling  effect  when  the 
spacing  between  array  elements  is  reduced.  We  let  the  radius 


of  the  circular  array  be  0.25  wavelength.  The  power  spectra 
when  there  are  two  incoming  signals  at  30  and  50°  are  plotted 
in  Figure  3.  Compared  with  the  previous  example,  we  see  that 
the  mutual  coupling  effect  is  much  stronger,  as  expected.  The 
uncompensated  curve,  which  ignores  the  mutual  coupling, 
shows  only  one  broad  peak.  The  CC  curve,  on  the  other 
f  hand,  still  correctly  resolves  the  two  DOAs. 

In  the  third  example,  we  examine  array  calibration  in  a 
very  noisy  environment.  We  reduce  the  signal-to-noise  ratio 
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<K°) 

Figure  3  MUSIC  power  spectrum  before  and  after  compensation  when  the  array  spacing  is  halved.  — Uncompensated, . CC 


to  0  dB.  Data  used  for  both  Atrue  and  C  contain  white  noise.  ment,  the  direction-finding  result  is  not  as  satisfactory.  The 

Two  signals  are  imposed  on  the  array,  and  the  MUSIC  power  uncompensated  result  is  again  the  poorest, 

spectra  are  plotted  in  Figure  4.  The  dashed  curve  is  the  In  the  fourth  example,  we  consider  the  antenna  gain  and 

uncompensated  result,  the  dash-dot  curve  is  the  CC  result,  delay  difference  between  the  array  elements,  in  addition  to 

and  the  solid  curve  is  the  reference  result.  Since  Atrue  is  the  mutual  coupling.  We  simulate  this  condition  in  NEC  by 
collected  in  the  presence  of  strong  noise,  the  reference  curve  using  different  load  impedances  for  different  antenna  ele- 

shows  the  degradation  of  the  MUSIC  algorithm  due  to  noise.  ments.  The  real  and  imaginary  parts  of  the  load  impedances 

The  CC  curve  includes  the  additional  degradation  from  esti-  are  randomly  selected  between  —  400  and  400  O.  The  power 

mating  the  coupling  matrix  using  noisy  data.  As  we  can  see,  spectra  for  one  signal  coming  at  0,  45,  or  90°  are  plotted  in 

when  the  coupling  matrix  is  estimated  in  the  noisy  environ-  Figure  5.  Some  of  the  uncompensated  results  show  false 


0  10  20  30  40  50  60  70  80  90 


,  0(°) 

Figure  4  MUSIC  power  spectrum  before  and  after  compensation  when  the  SNR  is  lowered  to  0  dB.  —Uncompensated,  . CC, 

— reference 
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Figure  5  MUSIC  power  spectrum  before  and  after  compensation  when  unknown  antenna  gains  and  delays  are  considered. 
—Uncompensated,  -• — CC 


peaks  in  the  spectrum.  The  peaks  corresponding  to  the  cor-  post  in  the  middle  of  the  array  to  simulate  a  real  antenna 

rect  DOAs  also  deviate  from  their  correct  locations.  This  array  layout.  The  post  is  0.2  wavelength  in  radius,  and  is 

implies  that  the  error  in  antenna  gain  and  delay  is  as  signifi-  approximated  by  36  vertical  wires  in  NEC.  We  again  consider 

cant  as  the  mutual  coupling  effect.  On  the  other  hand,  the  the  case  of  two  incoming  signals.  The  MUSIC  power  spectra 

sharp  peaks  from  the  CC  results  indicate  that  the  gain  and  are  plotted  in  Figure  6,  where  the  dashed  curve  is  the 

delay  errors  can  still  be  correctly  accounted  for  by  the  cou-  uncompensated  result,  the  dash-dot  curve  is  the  CC  result, 

pling  matrix  model.  This  is  not  surprising  since  the  condition  and  the  solid  curve  is  the  reference  result.  It  can  be  observed 

for  (1)  to  be  valid  is  still  satisfied.  that,  although  the  CC  result  is  much  better  than  the  uncom- 

Finally,  we  consider  an  example  where  a  platform  struc-  pensated  one,  it  cannot  achieve  the  same  resolution  as  the 

ture  is  placed  close  to  the  array.  We  put  a  center  conducting  reference  curve.  This  is  because  the  actual  array  response  is 


,<K°) 


Figure  6  MUSIC  power  spectrum  before  and  after  compensation  when  the  platform  effect  is  considered.  — Uncompensated, - CC, 

— reference 
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not  correctly  modeled  by  (1),  as  the  ideal  element  patterns 
used  in  (1)  do  not  include  the  effect  of  the  platform.  Thus, 
the  coupling  matrix  is  no  longer  angular  independent.  How¬ 
ever,  if  the  two  signals  are  not  too  closely  spaced  in  angle,  the 
DOA  estimates  are  still  satisfactory. 

IV.  CONCLUSIONS 

In  this  paper,  the  effect  of  mutual  coupling  on  direction-of- 
arrival  estimates  in  a  circular  array  is  simulated  using  rigor¬ 
ous  electromagnetic  computation.  The  compensation  of  the 
coupling  effect  using  the  coupling  matrix  approach  is  also 
examined.  Numerical  results  show  that  the  mutual  coupling 
effect  can  lead  to  significant  errors  in  direction  finding  when 
not  properly  accounted  for.  The  compensation  of  the  mutual 
coupling  effect  using  the  coupling  matrix  approach  is  found 
to  be  quite  satisfactory  in  most  cases.  The  coupling  matrix 
model  provides  a  very  sparse  description  of  the  array  re¬ 
sponse.  It  can  also  take  into  account  additional  antenna  gain 
and  delay  errors.  However,  the  performance  of  this  approach 
degrades  in  cases  when  the  array  calibration  data  contain  a 
high  noise  level  or  when  strong  platform  effects  are  present. 
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ABSTRACT:  A  frequency-interpolation  algorithm  for  electromagnetic 
scattering  data  from  large,  complex  targets  is  proposed.  It  is  based  on  the 
hybridization  of  the  existing  multiple-arrival  model  and  the  rational-func¬ 
tion  model.  At  each  point  on  the  target,  the  induced  current  is  parameter¬ 
ized  using  both  models  separately.  The  model  that  best  matches  the 
scattering  physics  at  that  location  is  chosen  based  on  the  resulting 
interpolation  error.  As  a  result,  those  regions  on  the  target  that  are  best 
described  by  ray-optical  phenomena  are  interpolated  by  the  multiple- 
arrival  model,  while  those  regions  that  exhibit  resonance  phenomena  are 
interpolated  by  the  rational-function  model.  Numerical  results  show  that 
the  hybrid  scheme  is  superior  to  either  the  multiple-arrival  model  or  the 
rational-function  model  for  a  target  containing  complex  features.  ©  2000 
John  Wiley  &  Sons,  Inc.  Microwave  Opt  Technol  Lett  27:  307-312, 
2000. 

Key  words:  frequency  interpolation;  method  of  moments;  multiple-arrival 
model;  rational-function  model 

1.  INTRODUCTION 

The  electromagnetic  scattering  data  from  a  complex  target 
are  often  of  interest  over  a  broad  band  of  frequencies.  To 
simulate  multiple-frequency  data  using  a  frequency-domain 
computational  electromagnetics  solver  is  an  exhaustive  proce¬ 
dure  as  the  computation  must  be  carried  out  one  frequency 
at  a  time.  The  problem  is  further  compounded  by  the  fact 
that  the  required  sampling  density  in  frequency  is  approxi¬ 
mately  proportional  to  the  size  of  the  target.  For  large 
complex  targets,  it  is  therefore  desirable  to  develop  algo¬ 
rithms  to  generate  a  dense  set  of  frequency  data  from  a  very 
sparse  set  of  computed  frequency  points  to  save  computation 
time.  This  frequency-interpolation  problem  was  recently  ad¬ 
dressed  in  [1]  using  a  model-based  approach  based  on  a 
high-frequency  model.  The  current  at  each  point  on  the 
target  is  modeled  by  a  collection  of  tinic-of-arrival  mccha- 


Contract  grant  sponsor:  Office  of  Naval  Research 
Contract  grant  number:  N000 14-98- 1-0178 

Contract  grant  sponsor:  Air  Force  MURI  Center  ior  Computational 
Electromagnetics 

Contract  grant  number:  AFOSR  F49620-96- 1-0025 


MICROWAVE  AND  OPTICAL  TECHNOLOGY  LETTERS  /  Vol.  27,  No  5,  December  5  2000 


nisms  arising  from  multiple  scattering  (see  Fig.  1).  The  time 
of  arrival  and  the  strength  of  the  different  mechanisms  are 
extracted  from  the  available  frequency  data  using  the  adap¬ 
tive  feature  extraction  (AFE)  algorithm  [2].  It  was  shown  that 
good  interpolation  results  can  be  obtained  even  for  very 
sparsely  sampled  data  in  frequency. 

In  this  paper,  we  set  out  to  improve  the  accuracy  of  the 
above  interpolation  algorithm  for  complex  targets  by  incorpo¬ 
rating  additional  scattering  physics  in  the  model.  The  multi¬ 
ple-arrival  model  is  well  matched  to  the  ray-optical  descrip¬ 
tion  of  fields  and  currents.  However,  more  general  targets 
contain  not  only  large-scale  features,  but  also  small  resonant 
features.  Although  these  resonant  contributions  are,  in  gen¬ 
eral,  weaker  than  the  returns  from  large-scale  features,  they 
do  impact  the  accuracy  of  the  interpolation  if  not  properly 
taken  into  account.  The  multiple-arrival  model,  however,  is 
not  a  good  model  for  resonance  phenomena  as  many  time- 
of-arrival  terms  are  required  to  adequately  describe  reso¬ 
nance.  A  much  better  model  is  the  rational-function  model, 
which  has  been  widely  used  in  frequency-interpolation  prob¬ 
lems  in  the  resonant  region  [3-8].  It  can  be  shown  that  the 
multiple-arrival  model  and  the  rational-function  model  are, 
in  fact,  complementary  models.  The  former  is  a  sum  of 
exponentials  in  frequency,  while  the  latter  is  a  sum  of  expo¬ 
nentials  in  time.  To  take  advantage  of  the  strength  of  each 
model,  we  propose  here  a  hybrid  interpolation  scheme  that 
uses  the  most  appropriate  model  over  different  portions  of 
the  target.  At  each  point  on  the  target,  we  parameterize  the 
induced  current  using  both  models  separately.  Next,  by  using 
the  resulting  parameterization  error  as  the  selection  crite¬ 
rion,  we  choose  the  model  that  best  matches  the  scattering 
physics  at  that  location.  The  total  scattered  field  is  then 
constructed  from  the  interpolated  currents.  Section  2  de¬ 
scribes  our  formulation.  The  two  models  and  the  methods  for 
parameterization  are  first  summarized.  They  are  followed  by 
the  hybrid  procedure.  Our  test  results  in  Section  3  show  that 
the  hybrid  scheme  is  superior  to  either  the  multiple-arrival  or 
the  rational-function  model  alone  for  a  target  containing 
complex  features. 


2.  FORMULATION 

2.1.  Multiple-Arrival  Model  and  AFE  Interpolation  Algorithm. 
In  the  multiple-arrival  model,  we  assume  that  the  induced 
current  at  a  point  on  the  target  can  be  written  as  a  summa¬ 
tion  of  multiple-scattering  mechanisms  from  M  different 
incident  paths  (see  Fig.  1): 

/(/>  =  L-4,exp(  -i2-.fi,)  P  (1) 

where  /  is  the  frequency,  is  the  time  of  arrival  for  mecha¬ 
nism  /,  and  At  is  the  corresponding  excitation  amplitude. 
This  time-of-arrival  model  is  based  on  high-frequency  ray- 
optical  phenomena,  and  has  been  well  utilized  in  the  electro¬ 
magnetics  community  [9-11].  The  last  term  in  this  model 
incorporates  an  additional  frequency-dependent  factor  ai9 
which  was  not  used  in  [1].  It  is  consistent  with  high-frequency 
diffraction  theory  [10,  12,  13],  and  improves  the  accuracy  of 
the  model.  Lastly,  the  extra  factor  in  the  denominator  is 
included  for  normalization. 

To  determine  the  unknowns  in  the  model  based  on  the 
available  frequency  samples  of  the  current,  we  use  the  adap¬ 
tive  feature  extraction  algorithm  [2].  The  parameterization  is 
carried  out  in  an  iterative  manner  starting  from  the  strongest 
mechanism.  The  sampled  frequency  response  is  first  pro¬ 
jected  onto  the  complex  conjugate  of  the  model  bases  for  all 
possible  values  of  tj  and  a;.  We  then  select  the  basis  that 
gives  the  maximum  projection  value.  This  is  described  as 
follows: 


At 


max( 


/,(/),  exp(/2irfy) 


(2) 


where  the  inner  product  in  the  above  formula  is  defined  as 
<«(/),  M/)>  =  ^  £a(f,)b(fp).  (3) 

P=  * 


After  the  strongest  feature  with  (Ahth  a,)  is  captured,  that 
feature  is  subtracted  from  the  signal  to  generate  a  remainder 
signal: 


Ji+ ,(/)  =  J,(f)  -  A,  exp(  —jl-irftj)  .  =  .  (4) 

The  above  procedure  is  iterated  to  extract  the  parameters  for 
each  scattering  mechanism  until  the  remaining  signal  reaches 
a  sufficiently  small  level. 

A  key  concept  in  AFE  is  to  use  random  frequency  sam¬ 
pling  during  the  original  data  generation  to  avoid  the  ambigu¬ 
ity  in  selecting  the  strongest  feature  [1].  Therefore,  the  fre¬ 
quencies  at  which  the  electromagnetic  computations  are 
carried  out  should  not  be  evenly  spaced.  In  addition,  since 
the  strongest  feature  contains  interference  from  the  weaker 
features,  the  amplitude  of  each  feature  is  not  extracted 
perfectly  using  AFE.  We  have  found  that,  by  performing  a 
linear  least  square  optimization  for  the  amplitude  parameter 
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after  the  whole  iteration  procedure,  the  AFE  performance  is 
significantly  improved. 


2.2.  Rational  Function  Model  and  Cauchy  Method.  The  multi¬ 
ple-arrival  model  and  AFE  can  be  used  to  effectively  carry 
out  frequency  interpolation  for  targets  containing  ray-optical 
phenomena.  For  target  resonance,  however,  it  is  more  appro¬ 
priate  to  use  a  rational-function  model  to  describe  the  fre¬ 
quency  resonance  mechanism.  The  rational-function  model  is 
introduced  below,  and  the  Cauchy  method  for  determining 
the  model  parameters  is  summarized.  The  reader  is  referred 
to  [5,  6]  for  a  detailed  description  of  the  method. 

Consider  the  frequency  response  of  a  system  H(s).  From 
circuit  theory,  it  is  customary  to  describe  H(s)  by  the  ratio  of 
two  polynomials  y4(s)  and  B(s)  as  follows: 


H(s)  = 


A(s) 

B(s) 


E  oksk 


k  =  0 

Q 


E  bksk 


k  =  0 


(5) 


For  Q  >  P,  this  is  a  pole  model,  and  thus  can  be  used  to 
describe  the  strong  resonance  behavior  of  the  current  in  the 
frequency  domain  by  the  proper  choice  of  the  pole  locations. 

Given  the  value  of  H(s)  and  its  frequency  derivatives  at 
some  frequency  points  sn,  the  Cauchy  problem  is  stated  as 

given  Hj(sn ),  for  j  =  0, . . . ,  /,  n  =  1, . . . ,  N 
find  P  and  Q ,  {ak,  k  =  and  {bk,  k  =  0,. . . , Q}. 

where  HKsn )  represents  the  yth  frequency  derivative  of  H  at 
sn.  Equation  (5)  can  be  rearranged  into  a  matrix  form  as 
follows: 


[A 


(6) 


where 


A{n,k) 


k\ 

(k  -  j)\ 


sl„k~j)u(k  -  j) 


(7) 


and 


j  k\ 

Bin,  k)  =  £  Cj'tHU-'KsJ— . —s^uik  -  /)■  (8) 

;  =  o 

Here,  u(k)  is  0  for  k  <  0  and  1  otherwise,  C}  =  (/!/(/!(/  - 
/)!)),  a  =  [ axa2a 3  aP]T,  and  b  =  [b{b2b3  •••  bQ]T.  To  solve 
for  the  polynomial  coefficients  ak  and  bk  from  the  matrix 
equation  (5),  singular  value  decomposition  is  utilized.  The 
model  orders  P  and  Q  can  be  estimated  from  the  number  of 
significant  singular  values.  The  singular  vector  corresponding 
to  the  smallest  singular  value  is  chosen  to  be  the  solution  of 
(6). 

Several  comments  are  in  order.  First,  the  quantity  to  be 
interpolated  in  our  problem  is  the  frequency  response  of  the 
current.  Thus,  H(s)  =  /(/)  where  s  =  j2irf.  Second,  we  only 
use  frequency  points,  and  not  frequency  derivatives  as  the 
input  to  the  Cauchy  algorithm.  This  is  for  ease  of  hybridiza¬ 
tion  with  the  AFE  algorithm.  Furthermore,  for  an  iterative 
solution  to  electromagnetic  integral  equations,  the  computa¬ 
tional  complexity  to  generate  the  frequency  derivative  of  the 


current  is  as  high  as  that  for  a  new  frequency  point.  Third, 
instead  of  estimating  the  orders  of  the  polynomial  functions, 
we  choose  Q  =  P  +  1  and  P  +  Q  +  1  =  N  in  our  implemen¬ 
tation.  We  find  that  our  results  are  not  very  sensitive  to  the 
order  estimation.  Lastly,  we  note  that  the  frequency  points 
need  not  be  equally  spaced  in  carrying  out  the  Cauchy 
method.  This  gives  us  the  possibility  to  use  the  same  set  of 
sampled  data  for  both  the  Cauchy  and  AFE  interpolation 
schemes. 

23.  Hybrid  AFE-Cauchy  Algorithm.  We  are  now  armed  with 
two  interpolation  algorithms.  The  AFE  algorithm  works  well 
for  low-order  multiple-scattering  mechanisms  that  resemble 
high-frequency  ray  optics,  while  the  Cauchy  algorithm  is 
better  suited  for  describing  resonance  phenomena.  The  es¬ 
sential  idea  of  the  hybrid  algorithm  is  to  choose  the  most 
appropriate  model  for  the  current  at  each  point  on  the  target. 
Our  approach  is  as  follows.  At  each  point  on  the  target,  we 
first  parameterize  the  induced  current  using  both  the  AFE 
and  Cauchy  algorithms  separately,  with  the  same  set  of 
sampled  frequency  responses.  The  model  that  best  matches 
the  scattering  physics  at  that  location  is  then  chosen.  In  the 
actual  implementation,  the  selection  is  made  by  comparing 
the  interpolated  data  from  the  two  models  against  the  brute- 
force  computation  at  several  new  frequency  points.  The  model 
with  the  smaller  interpolation  error  is  chosen  for  that  loca¬ 
tion.  As  a  result,  the  surface  of  the  target  can  be  divided  into 
two  regions — the  AFE  region  and  the  Cauchy  region.  Over 
the  AFE  region,  the  currents  exhibit  high-frequency  scatter¬ 
ing  phenomenology.  Over  the  Cauchy  region,  the  currents 
exhibit  strong  resonance.  Once  the  parameterization  and 
model  determination  have  been  made,  the  interpolated  cur¬ 
rents  on  the  target  over  the  frequency  band  of  interest  can  be 
generated.  Therefore,  the  scattered  far  field  can  be  obtained 
by  integrating  the  induced  current  over  the  target  at  any 
frequency  of  interest. 

It  has  been  pointed  out  that  the  AFE  algorithm  requires 
as  its  input  the  currents  computed  at  a  set  of  nonuniformly 
sampled  frequencies.  In  fact,  the  more  random  the  sampling, 
the  better  the  performance  of  the  AFE.  It  has  also  been 
noted  that,  in  the  Cauchy  algorithm,  the  frequency  points 
need  not  be  uniformly  sampled.  However,  equally  spaced 
frequency  points  increase  the  chance  for  the  Cauchy  method 
to  accurately  capture  the  poles  in  the  frequency  domain. 
Consequently,  we  find  contradictory  requirements  for  the 
optimal  sampling  in  frequency  between  the  AFE  and  Cauchy 
algorithms.  AFE  requires  random  sampling,  while  Cauchy 
prefers  uniform  sampling.  Here,  we  choose  a  compromise 
solution,  i.e.,  semirandom  sampling.  To  choose  the  N  fre¬ 
quency  points  needed  for  interpolation,  we  first  equally  divide 
the  whole  frequency  band  into  N  uniformly  spaced  subbands. 
The  sampling  point  within  each  subband  is  chosen  randomly 
based  on  a  uniform  probability  distribution.  The  sampled 
data  can  then  be  used  in  both  algorithms. 

3.  NUMERICAL  EXAMPLE 

We  consider  a  2-D  conducting  target  as  shown  in  Figure  2. 
Although  this  plate-like  target  looks  simple,  its  scattering 
characteristics  involve  a  variety  of  scattering  mechanisms  [14], 
An  £-poIarized  plane  wave  is  incident  at  an  angle  of  60°  from 
the  right.  The  reference  backscattering  data  versus  frequency 
are  first  computed  using  the  method  of  moments  (MoM)  at 
100  frequency  points  from  0.1  to  10  GHz  in  0.1  GHz  steps. 
The  frequency  response  of  the  backscattered  field  is  plotte 
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Figure  2  Geometry  of  the  conducting  plate  with  cavity  and  fin 


as  the  solid  line  in  Figure  3.  By  inverse  Fourier  transforming 
the  frequency  data,  the  range  profile  of  the  target  can  be 
generated.  It  is  plotted  in  Figure  4  as  the  solid  curve  (a 
Hamming  window  has  been  applied  to  reduce  range  side- 
lobes).  We  recognize  four  large  peaks  in  the  range  profile, 
which  correspond,  respectively,  to  the  scattering  from  the 
right  edge  of  the  plate,  the  aperture  of  the  partially  open 
cavity,  the  interior  of  the  cavity,  and  the  small  fin  at  the  left 
edge.  The  resonant  behavior  of  the  cavity  interior  can  clearly 


Figure  3  (a)  Comparison  of  the  backscattered  field  versus  fre¬ 
quency  between  the  brute-force  MoM  and  the  AFE  interpolation 
results  based  on  20  points,  (b)  Comparison  of  the  backscattered  field 
versus  frequency  between  the  brute-force  MoM  and  the  Cauchy 
interpolation  results  based  on  20  points,  (c)  Comparison  of  the 
backscattered  field  versus  frequency  between  the  brute-force  MoM 
and  the  hybrid  interpolation  results  based  on  20  points 


Figure  4  (a)  Comparison  of  the  range  profiles  generated  from  the 
brute-force  MoM  and  the  AFE-interpolated  data,  (b)  Comparison  of 
the  range  profiles  generated  from  the  brute-force  MoM  and  the 
Cauchy-interpolated  data,  (c)  Comparison  of  the  range  profiles  gen¬ 
erated  from  the  brute-force  MoM  and  the  hybrid-interpolated  data 


be  seen  as  dispersive  ringing  occurs  after  the  third  peak  and 
overlaps  with  the  fourth  peak.  We  should  note  that  the 
frequency  response  is  marginally  sampled  at  a  step  size  of  0.1 
GHz.  Further  undersampling  will  result  in  a  loss  of  features 
in  the  frequency  domain  or  aliasing  in  the  time/range 
domain. 

We  next  carry  out  the  frequency  interpolation  using  only 
20  frequency  points.  The  locations  of  the  20  frequency  points 
are  chosen  semirandomly,  as  discussed  in  the  last  section,  and 
are  indicated  as  small  arrows  along  the  frequency  axis  in 
Figure  3.  First,  the  target  current  is  interpolated  using  the 
AFE  algorithm  based  on  the  multiple-arrival  model  over  the 
whole  target  surface.  The  interpolated  result  is  plotted  as 
the  dashed  line  in  Figure  3(a).  Compared  to  the  reference 
result  from  the  brute-force  MoM,  AFE  gives  rather  good 
performance,  and  can  capture  most  of  the  features  in  the 
scattered  field.  The  largest  error  in  the  AFE  result  occurs  at 
frequencies  around  4.5  and  9  GHz,  which  are  close  to  the 
first  two  resonant  frequencies  of  the  small  cavity.  This  obser¬ 
vation  is  further  corroborated  in  the  corresponding  range 
profile  from  the  AFE-interpolated  data,  shown  as  the  dashed 
curve  in  Figure  4(a).  The  AFE  algorithm  performs  well  in 
predicting  the  three  strong  peaks  due  to  the  exterior  features 
on  the  target.  However,  it  has  trouble  predicting  the  resonant 
cavity  return,  which  is  rather  dispersive  in  range. 
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Next,  the  Cauchy  method  is  used  to  interpolate  the  cur¬ 
rent  over  the  entire  target  surface.  The  frequency  response 
and  the  range  profile  results  are  shown,  respectively,  in 
Figures  3(b)  and  4(b)  as  dashed  lines.  Since  the  rational 
function  is  a  pole  model,  we  find  that  the  Cauchy  algorithm 
tries  hard  to  capture  the  poles  of  the  system.  However,  since 
the  current  is  sparsely  sampled  in  frequency,  an  incorrect 
pole  location  and  pole  strength  are  generated.  This  is  particu¬ 
larly  true  for  low-order  multiple-scattering  events  that  have  a 
large  delay  spread  in  the  times  of  arrival.  The  long  delay 
spread  in  time  translates  into  rapid  oscillations  in  the  fre¬ 
quency  response,  and  a  very  high  model  order  is  needed  to 
adequately  model  the  response  using  poles.  From  the  range 
profile  in  Figure  4(b),  we  see  that  the  Cauchy  algorithm  gives 
extremely  poor  results.  Thus,  the  Cauchy  algorithm  cannot  be 
used  by  itself  to  carry  out  the  frequency  interpolation  of  large 
targets  that  are  dominated  by  ray-optical  phenomena. 

Finally,  we  apply  the  hybrid  interpolation  algorithm.  The 
selection  between  the  Cauchy-interpolated  and  the  AFE-in- 
terpolated  currents  at  each  point  on  the  target  surface  is 
made  by  comparing  them  against  the  reference  current  at 
three  extra  test  frequencies.  The  chosen  interpolated  cur¬ 
rents  are  then  used  to  generate  the  backscattered  field,  as 
shown  in  Figure  3(c).  When  compared  to  Figure  3(a)  and  (b), 
the  hybrid  result  shows  an  obvious  improvement  in  its  agree¬ 
ment  with  the  reference  result.  Note  that  the  error  near  the 
cavity  resonant  frequencies  in  Figure  3(a)  is  greatly  reduced 
in  Figure  3(c).  The  range  profile  in  Figure  4(c)  further  sub¬ 
stantiates  this  claim  as  both  the  scattering  centers  and  the 
resonance  of  the  cavity  are  well  predicted. 

Figure  5  shows  how  the  whole  target  is  automatically 
divided  into  two  regions  during  the  hybrid  procedure,  de¬ 
pending  on  the  local  scattering  physics.  The  Cauchy  region  is 
shown  within  the  open  lines,  while  the  AFE  region  is  shown 
within  the  solid  lines.  The  right  half  of  the  plate  fits  the 
multiple-arrival  model  very  well,  and  thus  is  chosen  as  the 
AFE  region.  The  cavity  in  the  middle  exhibits  strong  reso¬ 
nance,  and  can  be  well  parameterized  by  the  Cauchy  method. 
As  for  the  left  half  of  the  plate,  we  expect  it  to  be  an  AFE 
region  with  the  same  scattering  mechanisms  as  the  right-half 
part.  However,  since  the  delay  spread  in  the  times  of  arrival 
between  the  direct  incident  wave  and  the  scattered  wave  from 
the  fin  at  the  left  edge  is  small,  the  Cauchy  method  has  a 


Cauchy  region 


Figure  5  Illustration  of  how  different  portions  of  the  target  are 
parameterized  in  the  hybrid  approach.  TV  Cauchy  region  is  shown 
within  the  open  lines,  and  the  AFE  region  is  shown  within  the  solid 
lines 


sufficient  model  order  to  parameterize  the  scattering  behav¬ 
ior.  In  fact,  we  find  that  either  the  Cauchy  or  the  AFE  model 
may  be  chosen  in  this  region  with  equally  good  results. 

4.  SUMMARY 

In  this  paper,  we  have  proposed  a  frequency-interpolation 
algorithm  for  electromagnetic  scattering  data  from  large, 
complex  targets.  It  is  based  on  the  hybridization  of  the 
existing  multiple-arrival  model  and  the  rational-function 
model.  The  adaptive  feature  extraction  algorithm  is  used  to 
extract  the  parameters  in  the  multiple-arrival  model,  while 
the  parameters  in  the  rational-function  model  are  obtained 
using  the  Cauchy  method.  At  each  point  on  the  target,  we 
first  parameterize  the  induced  current  using  both  the  AFE 
and  Cauchy  algorithms  separately.  The  model  that  best 
matches  the  scattering  physics  at  that  location  is  automati¬ 
cally  chosen  based  on  the  resulting  interpolation  error.  Con¬ 
sequently,  those  regions  on  the  target  that  are  best  described 
by  ray-optical  phenomena  are  interpolated  by  the  multiple- 
arrival  model,  while  those  regions  that  exhibit  resonance 
phenomena  are  interpolated  by  the  rational-function  model. 
Numerical  results  for  a  target  containing  complex  features 
have  demonstrated  that  the  hybrid  model  results  in  more 
accurate  interpolation  than  either  of  the  two  models  alone. 
The  hybrid  interpolation  scheme  is  quite  robust,  and  can  lead 
to  significant  computational  savings  since  broadband  signa¬ 
ture  data  can  be  generated  from  a  very  sparse  set  of  com¬ 
puted  frequency  points.  Extension  of  this  algorithm  to  3-D, 
nonconducting  targets  is  currently  being  investigated. 
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Paper  [22] 

The  Effect  of  Mutual  Coupling  on  Direction  Finding  in  Smart  Antenna 

Applications 

Kapil  R.  Dandekar,  Hao  Ling  and  Guanghan  Xu 

Electromagnetic  computation  is  used  to  account  for  mutual  coupling  effects  in  a  smart 
antenna  system  and  improve  the  performance  of  direction  of  arrival  (DO A)  estimation. 
Experimental  results  of  the  Bartlett  and  Multiple  Signal  Classification  (MUSIC)  DOA 
algorithms  for  mobile  users  transmitting  at  1.8  GHz  are  given. 

Keywords:  Antenna  Arrays,  Array  Signal  Processing,  Mutual  Coupling 

Introduction:  The  effect  of  mutual  coupling  on  antenna  arrays  has  long  been  recognized 
in  the  array  processing  community  [1-5].  In  particular,  our  past  experimental  studies 
using  a  smart  antenna  testbed  showed  anomalous  DOA  results  [6].  For  a  single  mobile 
user  in  an  open  field  environment,  the  spatial  spectra  using  a  uniform  circular  array 
(UCA)  showed  a  large  sidelobe  not  corresponding  to  any  signal  multipath  or  interference. 

The  presence  of  this  sidelobe  offset  by  1 80°  from  the  direct  path  signal  component  led  to 
use  of  the  term  “rebounce”  sidelobe  to  describe  this  phenomenon.  Initial  speculation 
attributed  the  cause  of  rebounce  to  be  either  mutual  coupling  or  amplitude  and  phase 
mismatch  between  array  elements. 

In  this  letter,  we  set  out  to  explain  and  remove  the  rebounce  anomaly.  Our  approach  is  to 
compute  the  array  response  including  mutual  coupling  effects  using  the  computational 


electromagnetics  solver  Numerical  Electromagnetics  Code  (NEC)  [7].  This  computed 
array  response  is  then  applied  to  DOA  estimation  using  measurement  data  collected  from 
our  smart  antenna  testbed.  Our  results  show  that  the  rebounce  sidelobe  is  reduced 
significantly  when  mutual  coupling  is  taken  into  consideration. 

Array  Mathematical  Model:  Assume  that  there  are  M  elements  in  the  base  station 
antenna  array  located  at  (xj,  yO,  1  <  i  <  M.  A  steering  vector  characterizes  the  relative 
phase  response  of  each  antenna  array  element  to  an  incident  signal  with  DOA  0.  An 
expression  for  the  ith  component  of  the  theoretical  (“uncompensated”)  steering  vector  is 
given  by  Equation  (1): 

ai(0)  =  exp(-jk(xi  cos0  +  yj  sin0))  (1) 

where  k  is  the  wavenumber  of  the  incident  electromagnetic  radiation.  An  array  manifold 
matrix  contains  as  its  columns  a  collection  of  steering  vectors  corresponding  to  a  set  of 
angles.  DOA  algorithms,  such  as  Bartlett  and  MUSIC,  use  knowledge  of  the  steering 
vector  to  generate  a  spatial  spectrum  that  gives  signal  power  versus  direction.  The 
equations  for  the  Bartlett  and  MUSIC  spatial  spectra  are  given  respectively  as  [8]: 


^bartlett  (®) —  (*  (0)R5(0))/ (a  (0)a(9)) 

(2) 

JW  (e) = (SH  (e)a(e))/(5 "  (e)ni5(e)) 

(3) 

where  R  is  the  estimated  array  data  covariance  matrix  and  FT1  is  the  estimated  noise 

subspace  obtainable  through  an  eigenvalue  decomposition  of  R .  The  standard  approach 
to  DOA  estimation  is  to  use  the  ideal  steering  vector  given  in  Equation  (1)  in  the  above 
expressions.  However,  possible  interactions  among  the  antenna  elements  are  neglected. 


Array  Response  Including  Mutual  Coupling:  To  properly  account  for  mutual 
coupling,  NEC  simulation  is  used  to  compute  the  actual  array  response  (“compensated” 
steering  vectors)  to  incident  signals  from  each  direction  about  the  array.  The  array 
elements  are  each  four  collinear  half-wave  dipoles.  Each  of  the  dipole  elements  is 
represented  as  a  wire  divided  into  segments  for  NEC  analysis.  Lump  loading  is  used  to 
model  the  isolation  between  the  dipole  antennas  and  the  load  impedance  of  the  array 
elements.  To  simulate  the  array  response  to  a  signal  from  a  particular  DO  A,  the  induced 
currents  on  the  dipole  elements  due  to  a  plane  wave  incident  from  that  DOA  are  first 
computed.  The  array  response  is  the  vector  of  received  voltages  across  the  load 
impedance  of  each  antenna  element.  The  complete  array  manifold  containing  the  array 
response  from  all  simulated  incident  DOAs  is  used  when  performing  uplink  DOA 
analysis  with  measurement  data. 

Experiment  Setup:  Measurements  using  the  smart  antenna  testbed  were  made  in  an 
open  field.  A  7-element  UCA  with  an  operating  frequency  of  1 .88  GHz  was  used.  The 
immediate  environment  was  grassy  with  buildings  off  in  the  distance.  There  were  up  to 
two  signal  generators  with  dipoles  used  to  represent  mobile  users  transmitting  to  the 
basestation.  Data  was  collected  using  the  testbed  and  uplink  spatial  spectra  were 
considered  for  each  generator  transmitting  individually  and  both  generators  transmitting 
together. 

Results:  Figure  1  shows  the  spatial  spectrum  generated  using  the  Bartlett  DOA 
algorithm  for  measurement  data  due  to  a  single  mobile  user  located  at  48°.  The  dashed 


curve  is  generated  using  uncompensated  steering  vectors  and  shows  a  large  rebounce 
sidelobe  at  210°.  The  solid  curve  is  generated  using  the  NEC-computed  array  response 
that  takes  into  account  mutual  coupling.  While  the  effect  on  the  detected  DOA  and  main 
lobe  is  minimal,  the  rebounce  lobe  is  significantly  reduced  through  use  of  compensation. 
Figure  2  shows  the  MUSIC  spatial  spectrum  using  uncompensated  and  compensated 
steering  vectors  for  the  same  mobile  user.  Again,  the  rebounce  lobe  near  210°  is  reduced. 
In  addition,  the  main  lobe  indicating  the  actual  DOA  is  sharpened.  Measurements  at  a 
number  of  different  angles  about  the  array  nearly  always  verified  these  results  for  both 
DOA  algorithms. 

Figure  3  shows  the  uncompensated  and  compensated  MUSIC  spatial  spectra  for  a 
situation  in  which  there  are  two  mobile  users  (located  at  30°  and  78°)  transmitting 
simultaneously.  In  particular,  compensation  allows  the  DOAs  of  the  users  to  be  made 
much  more  prominent  compared  to  the  sidelobes.  For  all  tested  cases,  mutual  coupling 
compensation  typically  reduced  sidelobes  not  corresponding  to  mobile  users  by 
approximately  3  dB.  These  sidelobes  could  be  mistaken  for  either  another  mobile  user  or 
multipath  signal  energy. 

Conclusions:  The  anomalous  direction  finding  results  observed  in  a  smart  antenna 
testbed  has  been  explained  in  terms  of  mutual  coupling  effects.  This  work  shows  how 
off-line  calculations  using  the  widely  available  NEC  can  improve  direction  finding 
performance  in  smart  antenna  systems  without  the  need  for  costly  field  calibration.  The 


improvement  comes  in  the  form  of  reduced  rebounce  sidelobes  that  could  be  mistaken  for 
interference  or  multipath  signal  energy. 
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Figure  Captions 


Figure  1  -  Bartlett  Spatial  Spectrum  for  Single  Mobile  User  at  48° 

-  Compensated 

- Uncompensated 


Figure  2  -  MUSIC  Spatial  Spectrum  for  Single  Mobile  User  at  48° 

-  Compensated 

- Uncompensated 


Figure  3  -  MUSIC  Spatial  Spectrum  for  Two  Mobile  Users  at  30°  &  78° 

-  Compensated 

- Uncompensated 
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Tao  Su  and  Hao  Ling 

Dept,  of  Electrical  and  Computer  Engineering 
The  Univ.  of  Texas  at  Austin 
Austin,  TX  78712-1084 


Abstract 

The  coupling  matrix  is  a  standard  way  to  represent  the  mutual  coupling  effect 
between  elements  in  an  antenna  array.  In  this  paper,  two  commonly  used  approaches  to 
determine  the  coupling  matrix  are  examined  and  their  limitations  are  discussed.  The 
coupling  matrix  obtained  from  the  minimum  mean  square  error  matching  of  the  active 
and  stand-alone  element  patterns  is  shown  to  be  more  effective  than  the  mutual 
impedance  approach.  Furthermore,  the  coupling  matrix  model  is  extended  for  more 
complex  antenna  arrays  and  an  example  is  provided  to  illustrate  the  effectiveness  of  the 
extended  model. 


Key  Words:  array  mutual  coupling,  coupling  matrix 


I.  Introduction 


It  is  well  known  that  the  radiation/receiving  pattern  of  an  antenna  element  situated  in 
an  array  environment  can  be  quite  different  from  its  stand-alone  pattern  due  to  mutual 
coupling  effects.  The  antenna  element  pattern  in  the  array  environment  is  called  the 
active  element  pattern  in  the  antenna  community  [1].  In  applications  such  as  direction 
finding,  the  active  element  patterns  of  the  array  elements  are  needed  to  carry  out  direction 
of  arrival  estimation.  If  the  stand-alone  element  patterns  are  used  instead  of  the  active 
element  patterns,  significant  degradation  can  result  in  array  performance  [2]. 

To  properly  account  for  the  mutual  coupling  effect  and  for  simplicity,  it  is  customary 
to  express  the  active  element  patterns  as  the  product  of  a  coupling  matrix  and  the  stand¬ 
alone  element  patterns.  Although  the  coupling  matrix  concept  is  widely  utilized  for  array 
signal  processing,  two  very  different  approaches  have  been  adopted  in  the  literature  to 
interpret  and  determine  the  coupling  matrix.  The  first  approach  comes  out  of  the  antenna 
community.  Gupta  and  Ksienski  [3]  first  derived  the  coupling  matrix  from  a  microwave 
network  point  of  view.  The  array  is  viewed  as  a  multi-port  network  and  it  is  shown  that 
the  coupling  matrix  can  be  simply  related  to  the  generalized  impedance  matrix  of  the 
multi-port  network.  Since  the  generalized  impedance  matrix  is  an  intrinsic  property  of  the 
antenna  structure,  such  a  derivation  of  the  coupling  matrix  is  quite  attractive  from 
fundamental  electromagnetics  point  of  view.  This  approach  has  been  utilized  and 
extended  in  a  number  of  antenna  mutual  coupling  studies  [4-6]. 

The  second  approach  originated  from  the  array  signal  processing  community. 
Friedlander  and  Weiss  [7]  first  proposed  a  method  for  direction  finding  based  on  the 
coupling  matrix  model.  In  their  view,  the  coupling  matrix  is  simply  an  averaged  effect  of 


the  angular-dependent  relationship  between  the  active  element  patterns  and  the  stand¬ 
alone  element  patterns.  Based  on  this  assumption,  various  array  calibration  procedures 
have  been  developed  in  which  the  coupling  matrix  is  determined  by  minimum  mean 
square  error  (MMSE)  matching  of  the  two  sets  of  patterns  at  a  few  known  incident  angles 
[8,9]. 

While  both  of  these  approaches  have  been  used  with  success  on  arrays  consisting  of 
simple  dipole  elements,  the  assumptions  used  in  their  derivation  and  the  limitations  of 
each  approach  have  not  been  examined  in  detail  in  the  literature.  In  this  paper,  we  set  out 
to  compare  the  two  different  approaches  for  determining  the  coupling  matrix.  We  first 
review  the  two  approaches  and  discuss  the  limitations  of  each  approach  in  Section  II.  In 
Section  III,  numerical  results  are  generated  using  the  standard  numerical  electromagnetic 
solver  NEC  [9]  to  verify  our  observations.  In  Section  IV,  we  show  that  the  coupling 
matrix  model  can  be  further  extended  for  more  complex  antenna  arrays  and  an  example  is 
provided  to  illustrate  the  effectiveness  of  the  extended  model. 

II.  Mutual  coupling  models 

Assume  an  array  with  N  elements  is  located  in  the  xy-plane  and  receives  a  plane  wave 
incident  from  angle  $  in  the  azimuth  plane.  Let  us  consider  a  matrix  relationship  between 
the  active  element  patterns  and  the  stand-alone  element  patterns  as: 

A„W  =  MA<toW  (1) 

where  Atrue  and  A ,heo  are  the  active  and  stand-alone  element  patterns,  respectively,  and  M 
is  the  mutual  coupling  matrix.  A ,heo  has  N  rows  with  each  row  containing  the  stand-alone 


element  pattern  for  a  particular  element  in  the  array.  For  the  zth  row,  it  can  be  represented 
as: 

(2) 

where  f((j))  is  the  stand-alone  element  response  as  a  function  of  <j>  and  the  exponential 
term  is  the  explicit  position-dependent  phase  delay  for  the  Ith  element  located  at  (x„  y,). 
The  stand-alone  pattern  matrix  of  the  array  is  then  defined  as: 

A'heo  (0)  =  ta  (0)  a2  (<!>)■•■  aN  (<z>)]T  (3) 

The  active  pattern  matrix  Alrue  has  the  same  structure  as  A theo,  except  that  a,( <f>)  is  replaced 
by  the  active  element  pattern  of  the  zth  element.  The  two  different  approaches  in 
determining  the  mutual  coupling  matrix  M  are  discussed  below. 

A.  The  Z  approach 

In  the  first  approach,  the  antenna  array  with  N  elements  is  viewed  as  an  N  port 
microwave  network  [3].  The  relationship  between  the  antenna  terminal  voltages  and 
currents  on  reception  under  a  plane  wave  excitation  can  be  written  as: 

v\  =  Z\J\  +  •••  +  Z\ Ji  +  •••  +  ZWIN  +  Vloc 

Vt  —  Zixlx  +  •••  +  z,.7(.  +  •••  +  Z,,./v  +  Vioc  (4) 

V,  =  z  n\I\  +  •••  +  ZNiI,  +  •••  +  ZNNIN  +  VNoc 

In  the  above  expression,  Vioc  is  the  received  voltage  at  the  zth  port  when  all  the  antenna 

terminals  are  open  circuited,  and  Zy  is  the  mutual  impedance  between  the  zth  and  the  /h 

port  defined  as  Z,,  —Vjl ,  .If  we  assume  that  each  antenna  is  loaded  by  load 

r  ‘j  '/  J\/t=0,k*j  J 


impedance  ZiL  at  its  tenninals  and  use  the  relationship  Vi  =  -ZiLIi ,  we  can  express  the 
terminal  voltages  in  the  following  matrix  form: 

1+Zn _  h. L  ...  h. £L 

Zil  .  '  ZyL 

...  1+ZiL  ...  Zj!L 

Z  7  7 

}L  .  ,lL  . 

hi  ...  hi  ...  i  \  ha 

7  7  7 

L  ** iL  ^NL 

or 

ZV  =  V0C  (5) 

Note  that  if  we  stack  columnwise  the  received  voltage  vectors  V  obtained  for  different 
plane  wave  incident  angles,  the  resulting  matrix  is  in  fact  the  active  pattern  matrix  Atrue  in 
equation  (1)  (apart  from  an  unimportant  scaling  factor).  Similarly,  if  we  stack  the  open 
circuit  voltage  vectors  \oc  and  assume  that  they  are  the  same  as  the  voltages  received 
under  stand-alone  conditions,  the  resulting  matrix  is  the  stand-alone  pattern  matrix  A ,/,eo. 
Therefore,  the  mutual  coupling  matrix  in  (1)  can  be  written  as: 

M  =  Z'1  (6) 

This  completes  the  derivation  of  M  under  the  Z  approach. 

Since  the  normalized  impedance  Z  is  intrinsic  to  the  structure  and  is  completely 
independent  of  angle,  M  is  also  angle  independent.  This  is  a  highly  desirable  conclusion. 
However,  it  hinges  on  the  assumption  that  the  received  voltage  at  the  terminals  of  an 
antenna  when  all  the  antenna  elements  are  open-circuited  is  the  same  as  the  received 
voltage  when  all  the  other  elements  are  absent.  From  the  electromagnetic  point  of  view, 
this  is  clearly  not  true  as  the  induced  currents  on  the  antenna  elements  will  not  be 


identically  zero  under  the  open  circuit  condition.  As  a  result,  the  total  voltage  developed 
at  the  antenna  terminals  is  not  only  due  to  the  incident  field,  but  also  fields  produced  by 
the  induced  currents  on  other  elements.  Therefore,  while  equation  (5)  is  always  true,  the 
coupling  matrix  M  given  by  (6)  does  not  rigorously  relate  the  active  and  stand-alone 
element  patterns.  Under  certain  situations,  however,  this  model  may  work  approximately. 
For  example,  if  we  consider  an  array  comprising  of  half-wave  dipole  elements,  the 
induced  currents  on  every  element  should  be  very  small  under  the  open  circuit  condition 
and  the  Z  approach  may  be  quite  adequate.  This  will  be  examined  numerically  in  Section 

m. 


B.  The  C  approach 

In  the  second  approach,  the  mutual  coupling  matrix  is  viewed  as  an  average  of  the 
angular  dependent  relationship  between  the  active  and  stand-alone  element  patterns  [7]. 
Under  this  assumption,  M  can  be  obtained  by  using  the  MMSE  matching  between  the 
two  patterns  at  a  number  of  incident  angles: 

M=C  =  A„Ai„(A,toAl)-1  (7) 

We  refer  to  the  above  formulation  as  the  C  approach. 

The  effectiveness  of  the  C  approach  will  now  be  examined  from  the  induced  current 
point  of  view  based  on  the  method  of  moments  (MoM).  In  MoM,  the  integral  equation  of 
the  induced  currents  on  all  the  antennas  in  the  array  can  be  written  into  a  matrix  form  as 
follows: 


LI  =  V 


(8) 


where  I  is  the  current,  V  is  the  incident  plane  wave  excitation  and  L  is  the  moment 
matrix.  Suppose  we  remove  all  but  the  ith  antenna  in  the  array,  we  can  find  the  stand¬ 
alone  current  distribution  on  the  ith  antenna  by  solving  a  smaller  moment  system 
described  by: 

L°/I?  =  V,  (9) 

N  such  systems  can  be  written  for  all  the  elements  in  the  array.  Since  the  right  hand  side 
is  determined  by  the  incident  field  only,  it  does  not  change  whether  the  antenna  is  stand¬ 
alone  or  in  the  array  environment.  Thus  we  rewrite  equation  (8)  as 

'1,1  [L*  0 

[l  :  = 

_ij  L°  L°* 

or 

LI  =  L#I°  (10) 

By  solving  (10),  the  actual  current  distribution  is  related  to  the  stand-alone  current  by 

I  =  L-‘L°I°  (11) 

Next  we  try  to  extract  the  currents  at  the  terminals  of  the  antenna  elements  from  the  I 
vector  and  stack  them  over  all  the  incident  angles  to  arrive  at  the  active  pattern  matrix 
Alrue.  Similarly,  we  try  to  extract  the  terminal  currents  from  the  1°  vector  and  stack  them 
to  form  the  stand-alone  pattern  matrix  A ,heo-  First,  we  sift  out  the  terminal  currents  from  I 
and  the  corresponding  rows  from  the  matrix  (Lf’L0)  in  (1 1)  to  form: 

<12> 

where  the  subscript  ‘ term ’  indicates  the  rows  corresponding  to  the  antenna  terminals. 
Next,  we  sift  out  the  terminal  currents  from  1°  by  writing: 


(13) 


I°  = 

.  </T 

i - 

’/° 
x  I  .term 

i 

o 

s* 

1 - 

kn 

O 

I 

1 _ 

where  S,  is  the  shape  of  the  current  distribution  on  the  zth  element  and  I°Lterm  is  its 

corresponding  terminal  current.  Substituting  (13)  into  (12),  we  can  find  the  relationship 
between  the  actual  and  stand-alone  terminal  currents.  Note  that  the  moment  matrices  L 
and  L°  are  angular  independent.  If  the  shape  of  the  current  distribution  on  a  stand-alone 
antenna  element  is  also  angular  independent,  then  equation  (12)  can  be  rigorously  cast 
into  the  form  of  (1).  To  summarize,  we  have  shown  that  the  relationship  between  the 
active  element  pattern  and  the  stand-alone  element  pattern  is  angular  independent 
provided  that  the  shape  of  the  current  distribution  on  each  stand-alone  antenna  element  is 
independent  of  angle.  When  this  condition  is  satisfied,  the  mutual  coupling  matrix  can  be 
effectively  determined  by  (7). 

Let  us  now  examine  when  the  angular-independent  current  shape  condition  can  be 
met.  The  first  situation  has  been  discussed  in  [5],  where  all  the  antennas  are  vertical  wires 
and  the  incident  directions  have  the  same  elevation  angle.  In  this  case,  the  stand-alone 
current  distribution  of  an  antenna  is  exactly  the  same  over  all  the  incident  angles.  A 
second  common  situation  is  when  the  antenna  elements  are  working  near  resonance.  In 
this  case,  the  stand-alone  current  distribution  on  an  antenna  is  dominated  by  the  resonant 
mode  and  thus  has  approximately  the  same  shape  over  all  the  incident  angles.  An 
example  of  this  situation  is  an  array  of  half-wave  dipoles. 


III.  Numerical  examples 


In  this  section,  we  illustrate  the  conditions  for  the  Z  and  C  approaches  to  be  valid 
using  computer  simulation.  We  study  a  linear  array  consists  of  4  dipoles,  whose  centers 
are  spaced  0.45  wavelength  apart,  as  shown  in  Fig.  1.  The  simulations  are  carried  out 
using  NEC. 

In  the  first  example,  the  dipoles  are  half  wavelength  in  length.  All  the  dipoles  are 
loaded  with  73  ohms  of  load  impedance  at  the  center.  Elements  1  and  4  are  rotated  45° 
and  -45°,  respectively,  about  the  y  axis.  Elements  2  and  3  are  rotated  45°  and  -45°, 
respectively,  about  the  x  axis.  We  determine  the  coupling  matrix  using  both  the  C 
approach  and  the  Z  approach  for  this  structure.  The  current  at  the  center  of  each  antenna 
is  computed  for  incident  angles  from  -90°  to  90°  at  a  step  of  1°.  This  serves  as  the 
reference  active  element  pattern.  The  values  at  ±75°,  ±45°,  and  0°  are  used  to  calculate 
the  C  matrix  from  (7).  To  determine  the  Z  matrix,  we  first  obtain  the  mutual  admittance 
by  short  circuiting  all  the  antenna  ports,  driving  the  /h  antenna  with  a  voltage  source,  and 

computing  Yij  =  h/Vj\v  •  By  inverting  the  resulting  generalized  admittance  matrix, 

we  get  the  mutual  impedance  Zy.  The  normalized  impedance  matrix  Z  is  then  obtained 
from  (5). 

Once  we  obtain  the  two  coupling  matrices,  we  use  them  to  calculate  the  active 
element  patterns  from  the  stand-alone  element  patterns.  Notice  that  since  the  elements  are 
tilted,  even  the  stand-alone  element  patterns  are  not  isotropic.  The  active  element  patterns 
of  antenna  elements  1  and  2  are  plotted  in  Figs.  2(a)  and  2(b),  respectively.  The  solid 
curves  are  the  active  element  patterns  directly  computed  by  NEC.  The  dashed  curves  are 
obtained  using  the  Z  approach  and  the  dotted  lines  are  modeled  using  the  C  approach. 
We  find  that  the  dashed  curves  are  quite  close  to  the  computed  ones,  but  not  exactly  the 


same.  This  is  because  the  Voc  condition  described  in  Section  IIA  is  only  approximately 
satisfied,  even  for  an  array  of  half-wave  dipoles.  On  the  other  hand,  the  curves  using  the 
C  approach  match  almost  perfectly  with  the  computed  ones.  In  this  case,  the  antennas 
operate  in  resonance  and  the  stand-alone  current  distribution  has  the  same  shape  over  all 
incident  angles,  even  though  they  are  not  all  vertically  aligned. 

To  further  observe  the  effectiveness  of  the  mutual  coupling  models,  we  apply  the 
active  element  patterns  to  a  direction  finding  problem,  which  is  quite  sensitive  to  the 
element  patterns.  Two  uncorrelated  incoming  signals  are  assumed  from  the  incident 
angles  of  50°  and  65°.  The  signal  to  noise  ratio  (SNR)  is  set  to  30  dB.  The  angles  of 
arrival  are  determined  using  the  standard  MUSIC  algorithm,  and  the  resulting  normalized 
power  spectra  are  plotted  in  Fig.  3.  The  solid  curve  is  the  result  from  the  reference  active 
element  patterns.  For  comparison,  the  result  by  neglecting  mutual  coupling  effect  is 
obtained  by  using  the  stand-alone  element  patterns  and  is  plotted  as  the  dash-dot  curve. 
The  large  degradation  shows  the  importance  of  taking  mutual  coupling  into 
consideration.  The  dashed  curve  is  the  result  of  the  Z  approach.  The  dotted  curve  is 
obtained  using  the  C  approach.  It  is  observed  that  the  C  approach  performs  almost  as 
well  as  the  reference  active  element  patterns,  while  the  Z  approach  leads  to  some 
degradation  in  the  peak  position  and  width  of  the  power  spectra. 

In  the  second  example,  we  replace  the  half-wave  dipole  elements  with  full-wave 
dipoles.  Compared  with  the  previous  example,  the  Voc  condition  is  strongly  violated  for 
the  Z  approach  and  the  near-resonant  condition  is  not  satisfied  for  the  C  approach.  Fig.  4 
shows  the  computed  and  modeled  active  element  patterns.  This  time,  neither  the  Z  nor 
the  C  approach  is  close  to  the  reference  result. 


Finally  we  make  the  full  wave  dipoles  perpendicular  to  the  xy  plane.  Now  the  vertical 
wire  condition  is  satisfied  for  the  C  approach.  The  resulting  element  patterns  are  plotted 
in  Fig.  5.  The  Z  approach  still  works  poorly  as  in  example  2.  However,  the  C  approach 
perfectly  models  the  active  element  patterns.  As  discussed  in  the  Section  IIB,  the 
matching  in  this  case  is  exact. 

IV.  Extension  of  the  mutual  coupling  model 

As  can  be  seen  in  the  examples,  the  C  approach  is  more  accurate  than  the  Z  approach. 
Yet  the  situations  when  the  C  approach  can  be  applied  is  still  quite  limited.  The  approach 
fails  for  more  complex  antenna  structures.  By  observing  equation  (11)  we  can  see  that  if 
the  array  elements  consist  of  thin  wires,  each  piece  of  the  wires  affects  the  active  element 
pattern  in  the  same  manner  as  a  dipole  antenna.  Thus  for  certain  cases,  the  model  in  (1) 
can  be  extended  to  include  the  coupling  from  all  parts  of  the  array.  For  example,  if  we 
have  an  array  of  Yagi  antennas,  we  may  include  the  coupling  from  the  parasitic  elements 
to  model  the  active  element  patterns  more  accurately.  More  specifically,  we  extend  the 
model  as  follows: 


where  A \heo  is  the  stand-alone  receiving  patterns  of  the  parasitic  elements  and  C' 
describes  the  coupling  from  the  parasitic  elements  to  the  driven  elements.  Like  the 
standard  C  approach,  (14)  is  valid  if  all  the  parasitic  elements  are  vertical  or  their  lengths 
are  close  to  resonant  length.  The  MMSE  matching  can  again  be  used  to  determine  the 
coupling  matrix,  except  that  the  number  of  incident  angles  needed  to  give  a  unique 


solution  is  increased. 


As  an  example,  we  look  at  an  array  consisting  of  4  parallel  Yagi  antennas  as  shown  in 
Fig.  6(a).  All  the  elements  in  the  array  are  perpendicular  to  the  xy-plane  and  the 
separation  between  two  adjacent  Yagis  is  1  wavelength.  Each  Yagi  antenna  consists  of  3 
wires  with  the  center  wire  being  the  driven  element  as  shown  in  Fig.  6(b).  The  stand¬ 
alone  element  pattern  is  plotted  in  Fig.  6(c).  Next,  the  coupling  matrix  in  (14)  is 
determined  from  the  computed  active  element  patterns  at  12  incident  angles  uniformly 
distributed  between  -180°  and  180°.  The  reference  and  the  modeled  active  element 
patterns  of  the  first  and  second  antenna  are  plotted  in  Figs.  7(a)  and  7(b),  respectively. 
The  solid  curves  are  the  reference  patterns  and  the  dotted  curves  are  obtained  using  the 
extended-C  approach.  They  match  with  the  reference  patterns  very  well.  As  a 
comparison,  we  also  model  the  active  element  patterns  using  the  Z  approach  and  the  C 
approach.  The  results  are  shown  in  dashes  and  dash-dot  curves,  respectively.  The 
agreement  is  obviously  not  as  good.  This  is  because  the  open-circuit  voltage  condition 
does  not  hold  for  the  Z  approach  and  the  constant  current  shape  condition  for  each  array 
element  is  violated  for  the  C  approach. 

V.  Conclusions 

In  this  paper,  the  two  commonly  used  approaches  for  determining  the  mutual 
coupling  matrix  in  array  antennas  have  been  discussed.  It  has  been  shown  that  the  Z 
approach  is  an  approximation  that  applies  when  the  element  pattern  under  the  open 
circuit  condition  can  be  well  approximated  by  the  stand-alone  element  pattern.  It  has  also 
been  shown  that  the  C  approach  holds  true  whenever  the  current  distribution  of  a  stand¬ 
alone  element  has  the  same  shape  over  all  incident  angles.  This  condition  is  satisfied  for 


vertical  wire  antennas  or  near-resonant  antennas.  We  have  found  from  our  numerical 


examples  that  the  C  approach  is  in  general  more  accurate  than  the  Z  approach.  We  have 
also  extended  the  C  approach  to  model  more  complex  wire  elements  by  considering  each 
wire  as  an  individual  antenna.  An  example  of  an  Yagi  array  has  been  provided  to 
demonstrate  the  effectiveness  of  the  extended-C  model. 
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Fig.  1 .  Array  Geometry 
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Fig.  4.  Computed  and  modeled  active  element  patterns  for  array  of  rotated 
full  wave  dipoles. 
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Fig.  5.  Computed  and  modeled  active  element  patterns  for  array  of  vertical 
full  wave  dipoles. 
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Fig.  7.  Computed  and  modeled  active  element  patterns  for  an  array  of  Yagi 
antennas. 
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Abstract:  The  genetic  algorithm  (GA)  is  used  to  design  patch  shapes  for 
microstrip  antennas  on  FR-4  substrate  for  broadband  applications.  Measurement 
results  of  the  GA-optimized  designs  show  good  agreement  with  numerical 
prediction.  The  optimized  patch  design  achieves  a  four-fold  improvement  in 
bandwidth  when  compared  to  a  standard  square  microstrip  antenna. 

Introduction:  It  is  well  known  that  standard  microstrip  patch  antennas  exhibit 
very  narrow  bandwidth.  Various  broadbanding  methods  have  been  proposed  to 
date.  For  instance,  adding  parasitic  patches,  using  thick  air  substrate,  stacking 
patches  and  using  shorting  post  for  reactive  loading  are  well  known  techniques  to 
extend  microstrip  bandwidth  [1].  Recently,  Johnson  and  Rahmat-Samii  reported 
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on  the  use  of  the  genetic  algorithm  (GA)  'to  search  for  novel  patch  shapes  for 
broadband  operation  [2].  The  attractiveness  of  GA  shape  optimization  is  that 
improved  bandwidth  performance  can  be  achieved  without  increasing  overall 
volume  or  manufacturing  cost.  They  used  thick  air  substrate,  and  explored  metallic 
patch  sizes  up  to  half  a  wavelength. 

In  this  paper,  we  also  examine  the  use  of  GA  for  broadband  applications.  In 
contrast  to  the  work  of  Johnson  and  Rahmat-Samii,  we  employ  FR-4  as  the 
substrate  material,  since  it  is  the  most  commonly  used  material  in  wireless  devices. 
In  addition,  fewer  geometrical  constraints  are  used  in  the  GA  in  hope  of  obtaining 
better  global  optimum.  We  report  a  four-fold  bandwidth  improvement  from  our 
GA-optimized  microstrip  shape  as  compared  to  a  standard  square  microstrip 
antenna. 

GA  Optimization:  GA  is  implemented  to  optimize  the  microstrip  patch  shape  in 
order  to  achieve  broad  bandwidth.  In  our  GA,  we  use  a  two-dimensional  (2-D) 
chromosome  to  encode  each  patch  shape  into  a  binary  map  [3].  The  metallic  sub¬ 
patches  are  represented  by  ones  and  the  no-metal  areas  are  represented  by  zeros. 
Since  it  is  more  desirable  to  obtain  optimized  patch  shapes  that  are  well-connected 
from  the  manufacturing  point  of  view,  a  2-D  median  filter  is  applied  to  the 
chromosomes  to  create  a  more  realizable  population  at  each  generation  of  the  GA. 
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To  evaluate  the  performance  of  each  patch  shape,  a  full-wave 
electromagnetic  patch  code  is  used  to  predict  its  bandwidth  performance  [4].  The 
formulation  of  the  code  is  similar  to  that  described  in  [5]  and  is  based  on  the 
solution  to  the  electric  field  integral  equation  with  the  periodic,  layered  medium 
Green’s  function  as  its  kernel.  Roof-top  basis  functions  are  used  to  expand  the 
unknown  current  on  the  metal  patch  and  fast  Fourier  transform  is  used  to  accelerate 
the  computation  of  the  matrix  elements.  Because  of  the  assumed  periodicity  of  the 
patches  in  this  code,  we  use  a  large  enough  period  to  simulate  a  single  patch. 

The  design  goal  is  to  broaden  the  bandwidth  of  a  microstrip  antenna  with  a 
center  frequency  of  2  GHz  by  changing  the  patch  shape.  To  achieve  the  design 
goal,  the  cost  function  is  defined  as  the  average  of  those  Sn  values  that  exceeds 
-KMB  (i.e.,  VSWR=2:1)  within  the  frequency  band  of  interest.  The  target 
frequency  range  is  between  1 .9  GHz  and  2. 1  GHz. 

Based  on  the  cost  function,  the  next  generation  is  created  by  a  reproduction 
process  that  involves  crossover,  mutation  and  2-D  median  filtering.  A  2-point 
crossover  scheme  using  three  chromosomes  is  devised.  The  process  selects  three 
chromosomes  as  parents,  and  divides  each  chromosome  into  three  parts. 
Intermingling  the  three  parent  chromosomes  then  makes  three  child  chromosomes. 
This  crossover  scheme  exhibits  a  more  disruptive  characteristic  for  regeneration 
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than  the  conventional  1 -point  or  2-point  crossover.  It  serves  to  counteract  against 
the  median  filtering  effect  and  is  found  to  result  in  better  convergence  rate.  The 
reproduction  process  is  iterated  until  the  cost  function  converges  to  a  minimum 
value. 

Results:  Fig.  1(a)  shows  the  shape  of  the  GA-optimized  microstrip.  A  72mm  x 
72mm  square  design  area  in  which  the  metallic  patch  can  reside  is  discretized  into 
a  16x16  grid  for  the  chromosome  definition.  The  thickness  of  the  FR-4  substrate 
(dielectric  constant  of  4.3)  is  1 .6  mm.  In  the  GA-optimized  shape,  the  gray  pixels 
are  metal  and  the  white  pixels  have  no  metal.  The  white  dot  shows  the  position  of 
the  probe  feed.  To  experimentally  verify  the  GA  design,  we  have  built  and 
measured  such  a  microstrip  patch.  Fig.  1(b)  shows  the  return  loss  comparison 
between  the  measurement  and  simulation  results.  The  solid  line  is  the  measurement 
result  taken  from  an  HP 8 753 C  network  analyzer.  The  square  dots  represent  the 
simulation  result.  Good  agreement  can  be  observed  between  the  measurement  and 
simulation.  The  graph  shows  a  bandwidth  of  6.16%  by  simulation  and  6.18%  by 
measurement.  This  is  about  three  times  that  of  a  square  microstrip  antenna  (36  mm 
x  36  mm),  which  has  a  bandwidth  of  1.98%.  Further  improvements  in  the 
bandwidth  can  be  obtained  from  the  GA  by  increasing  the  grid  resolution  from 
16x16  to  32x32.  Figs.  2  (a)  and  2  (b)  show  respectively  the  GA-optimized  patch 
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shape  and  the  bandwidth  performance  in  the  higher  resolution  design.  The 
bandwidth  is  found  to  be  8.04%  by  simulation  and  8.10%  by  measurement.  This  is 
about  four  times  that  of  a  square  microstrip  antenna. 

Finally,  the  operating  principle  of  the  GA-optimized  shape  is  interpreted.  It 
is  clear  from  the  two  frequency  dips  in  Fig.  2(b)  that  the  antenna  contains  two 
operating  modes  that  are  very  closely  spaced  in  frequency.  We  have  verified  the 
two  modes  by  examining  the  current  distributions  on  the  patch  at  1.99  and  2.07 
GHz.  In  addition  to  the  dual-mode  principle,  another  important  bandwidth 
enhancement  effect  is  achieved  through  the  ragged  edge  shape.  We  have  found 
that  when  the  patch  is  restricted  to  single-mode  operation  (by  imposing  symmetry 
constraints),  the  introduction  of  ragged  edges  in  the  GA-optimized  shape  can 
enhance  the  bandwidth  by  about  30%.  Therefore,  the  GA-optimized  design 
combines  both  the  dual-mode  operation  and  ragged  edge  shape  to  achieve  the 
broadest  bandwidth. 

Conclusions:  Optimized  patch  shapes  for  microstrip  antennas  on  thin  FR-4 
substrate  have  been  investigated  using  the  genetic  algorithm.  The  optimized  shape 
shows  a  four-fold  improvement  in  bandwidth  when  compared  to  a  standard  square 
microstrip  antenna.  This  result  has  been  verified  by  laboratory  measurement.  The 
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basic  operating  principle  of  the  optimized  shape  can  be  explained  in  terms  of  a 
combination  of  dual-mode  operation  and  ragged  edge  shape. 
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16x16 

(a) 


CHI  A/R  loo  MAG  3  dB/  REF  0  dB  3;  -12. 169  dB 


(b) 

Fig.  1.  (a)  GA-optimized  microstrip  antenna  using  16x16  resolution  within  a 
72mm  x  72mm  area.  The  gray  pixels  are  metal  and  the  white  dot  shows  the 
position  of  the  probe  feed,  (b)  Return  loss  (dB)  of  the  GA-optimized  antenna  from 
simulation  (-)  and  measurement  (□). 
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(a) 


CHI  A/R  log  MAO  2.S  dB/  REF  0  OB  3;  -10.374  dB 


(b) 


Fig.  2.  (a)  GA-optimized  microstrip  antenna  using  32x32  resolution  within  a 
72mm  x  72mm  area.  The  gray  pixels  are  metal  and  the  white  dot  shows  the 
position  of  the  probe  feed,  (b)  Return  loss  (dB)  of  the  GA-optimized  antenna  from 
simulation  (-)  and  measurement  (□). 
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A  wavelet-based  method  is  proposed  to  effectively  precondition  3-D 
electromagnetic  integral  equations.  The  approximate-inverse 
preconditioner  is  constructed  in  the  wavelet  domain  where  both  the 
moment  matrix  and  its  inverse  exhibit  sparse,  multilevel  finger  structures. 

The  inversion  is  carried  out  as  a  Forbenius-norm  minimization  problem. 
Numerical  results  on  a  3-D  cavity  show  that  the  iteration  numbers  are 
significantly  reduced  with  the  preconditioned  system.  The  computational 
cost  of  the  preconditioner  is  kept  under  0(NlogN). 

Introduction:  There  is  growing  interest  in  the  computational  electromagnetics  community 
on  the  use  of  the  fast  multipole  method  (FMM)  [1]  for  solving  large-scale 
electromagnetic  integral  equations.  With  the  multi-level  implementation  of  the  FMM, 
the  computational  complexity  and  storage  requirement  can  be  reduced  to  OiNlogN)  for 
the  (moment  matrix)-(vector)  product.  In  order  to  take  advantage  of  this  reduced 
complexity,  iterative  solvers  must  be  used.  However,  the  convergence  rate  of  the 
iterative  solution  process  is  strongly  problem-dependent.  For  scattering  geometries 
involving  multiple  interactions  such  as  partially  open  cavities,  the  moment  matrices  are 
severely  ill-conditioned  and  preconditioning  is  need  to  accelerate  the  convergence  rate 
[2].  Some  well-documented  preconditioning  methods  such  as  incomplete  LU- 
factorization  (ILU)  [3]  may  be  effective,  but  usually  require  well-above  OiNlogN) 


complexity  to  implement.  In  this  work,  we  propose  a  new  wavelet-based  method  to 
construct  an  effective  preconditioner  for  3-D  moment  equations.  The  algorithm  is  based 
on  determining  an  approximate  inverse  for  the  moment  matrix  in  the  wavelet  transform 
domain.  The  computational  cost  of  the  preconditioner  is  kept  under  0(NlogN),  making  it 
compatible  with  the  FMM. 

Problem  Formulation:  Consider  a  3-D  moment  equation  in  the  space  domain: 

[Z]  J  =  E  (1) 

where  [Z],  J,  and  E  are  respectively  the  moment  matrix,  the  induced  current  vector,  and 
the  excitation  vector.  If  (1)  is  left-preconditioned  by  a  preconditioner  [P],  we  obtain: 

[P][Z]J  -  [P]E  (2) 

Our  approach  is  to  try  to  construct  a  [P]  that  is  the  approximate  inverse  of  [Z]  using  the 
wavelet  transform.  The  moment  equation  in  (1)  can  be  represented  using  wavelet  basis 

as:  [Z]  J  =  E  (3) 

where 

[Z]  =[M]t[Z][M],  [J]=[M]tJ,  [EJ  =[M]tE,  (4) 

and  [M]  is  the  unitary  wavelet  transform  matrix.  Note  that  in  3-D  problems  [Z]  is 
strongly  diagonal-dominant.  Consequently,  [Z]  can  be  effectively  approximated  by  a 
sparse  matrix  with  the  multi-level  “finger”  pattern  discussed  in  [4].  Furthermore,  we 
observe  from  (4)  that: 

[Zr,=[M]T[Z]-1[M]=[Zrl  (5) 

Since  the  smooth  parts  in  [Z]'1  are  converted  into  small  wavelet  coefficients  through  the 
transform,  [Z]  ‘can  again  be  approximated  by  a  sparse  matrix  with  the  same  finger 


pattern  as  [Z] .  With  the  chosen  sparsity  pattern  for  both  [Z]  and  [Z]'1 ,  the  approximate- 

inverse  preconditioner  for  [Z]  can  now  be  solved  efficiently  by  casting  it  as  a  Forbenius- 
norm  minimization  problem  [5]: 


^n||[p][Z]-/F 


where  /  is  the  identity  matrix.  Since 


where  p.  and  ej  are  respectively  the y'-th  rows  of  [P]  and  /,  the  solution  of  (6)  becomes 
the  following  N  independent  least  squares  problems: 

MmllPjfZJ-ej  \  j  =  1,2,..., N  (8) 


Once  [P]  is  solved  for  in  the  wavelet  domain,  the  approximate-inverse  preconditioner  [P] 
can  be  obtained  by  the  inverse  wavelet  transform: 

[P]  =  IM][PJ[M]t  « [Z]'1  (9) 

Combining  (2)  and  (9)  we  have  the  following  preconditioned  moment  equation: 

[M][P][M]t[Z]J  =[M][P][M]tE  (10) 

To  summarize,  the  wavelet  preconditioner  [P]  is  generated  by  first  finding  the 
approximate  wavelet  transform  of  the  moment  matrix  [Z]  using  (4),  and  then  solving  the 
Forbenius-norm  minimization  problem  in  (8).  Once  [P]  is  obtained,  the  preconditioning 
operation  is  carried  out  by  a  series  of  matrix-vector  multiplication  operations  in  (10).  We 
now  consider  the  computational  cost  of  constmcting  [P]  and  that  for  the  preconditioning 

operation.  To  construct  [P],  we  must  obtain  an  approximate  [Z].  Due  to  the  strong 


source  singularity  in  3-D  problems,  [Z]  can  be  made  approximately  sparse  with  0(N) 
non-zero  elements  by  ignoring  the  far-field  interactions.  Therefore,  the  computation  of 

[Z]  can  implemented  with  fast  2-channel  filterbank  filtering  with  about  O(N)  operations. 

Once  [Z]  is  obtained,  the  complexity  to  solve  (8)  can  be  made  proportional  to  the 

problem  size  N,  provided  that  we  carefully  control  the  sparsity  patterns  of  [Z]  and  its 
inverse.  Next  we  consider  the  complexity  to  carry  out  the  preconditioning  operation  on 
the  left-hand  side  of  (10).  Since  there  are  only  0(NlogN)  non-zero  elements  in  the 
wavelet  transform  matrix  [M],  the  product  of  [M]  with  any  vector  requires  0(NlogN) 
operations.  Note  also  that  [P]  is  a  sparse  matrix  with  0(N)  non-zero  elements. 
Therefore,  if  the  preconditioning  operation  is  carried  out  as  a  series  of  matrix-vector 
products  on  [Z]  J  or  E,  the  total  computation  cost  is  limited  to  0(NlogN)  per  iteration. 

Numerical  Results:  The  algorithm  is  tested  using  a  3-D  conducting  rectangular  cavity 
with  an  open  end.  The  (width)x(height)x(depth)  dimensions  of  the  cavity  are  RxRx\.5R 
where  R  is  a  size  parameter.  The  incident  plane  wave  is  horizontally  polarized  with  a 
frequency  of  5.9GHz  and  makes  an  angle  of  45°  from  the  cavity  opening  in  the  elevation 
plane.  The  problem  is  formulated  in  terms  of  an  electric-field  integral  equation.  The 
Rao-Wilton-Glisson  basis  functions  are  used  to  form  the  original  space-domain  moment 
equation  with  a  discretization  size  of  about  Ay  10  (or  0.2").  The  grayscale  magnitude 

plots  of  the  wavelet-based  moment  matrix  [Z]  and  its  inverse  [Z]'1  are  shown 
respectively  in  Figs.  2(a)  and  2(b)  with  a  problem  size  of  iV=1024  (/?~2.1").  The  wavelet 
filter  used  in  the  transform  is  the  Daubechies  filter  of  order  6,  and  the  maximum  wavelet 
transform  level  is  Z,=4.  As  expected,  both  matrices  exhibit  the  multi-level  “finger” 


structure.  Accordingly  we  propose  to  use  a  fixed  sparse  matrix  pattern  shown  in  Fig.  2 
with  finger  width  D  for  both  [Z]  and  [Z]'1  in  solving  the  Frobenius-norm  minimization 
problem.  If  the  parameters  D  and  L  are  kept  constant,  the  complexity  to  solve  (10)  is 
proportional  to  the  problem  size  N.  To  investigate  the  performance  of  the  preconditioner, 
we  proportionally  increase  the  physical  size  of  the  cavity  with  the  parameter  R  changing 
from  1.2"  to  5.1"  so  that  N  increases  from  256  to  8192,  while  the  discretization  interval 
and  the  frequency  remain  unchanged.  The  resulting  iteration  numbers  in  solving  the 
moment  equations  as  a  function  of  the  problem  size  are  shown  in  Fig.  3.  The  iterative 
solver  is  BICGSTAB  and  a  relative  residual  error  of  0.001  is  used  as  the  stopping 
criterion.  The  maximum  wavelet  transform  level  is  L=5  and  the  finger  width  D  of  the 
sparsity  pattern  is  16  for  N<1024  and  32  for  N>1024.  We  observe  that  with 
preconditioning,  the  iteration  numbers  are  very  small  and  grow  at  a  very  slow  rate  with 
respect  to  the  problem  size.  The  results  demonstrate  the  effectiveness  of  the  new 
wavelet-based  preconditioner  for  3-D  moment  equations. 

Conclusions :  We  have  proposed  a  wavelet-based  preconditioner  for  3-D  electromagnetic 
integral  equations  in  this  work.  Numerical  results  showed  the  preconditioner  is  very 
effective  for  ill-conditioned  cavity  structures.  The  total  computational  cost  for  the 
preconditioning  can  be  kept  under  0(NlogN).  The  new  algorithm  is  compatible  with  fast 
boundary-integral  algorithms  such  as  the  multilevel  FMM. 
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Figure  1 .  (a)  The  wavelet  transformed  moment  matrix  [Z]  for  the  rectangular  cavity  in 


logarithmic  scale  with  -£,=4(77=1024).  (b)  The  inverse  of  the  transformed  matrix  [Z]'1 . 


Figure  2.  The  non-zero  pattern 
used  to  solve  the  Frobenius- 
norm  minimization  problem. 


Figure  3.  The  iteration  numbers  vs.  problem  sizes  for 
the  cavity  problem.  ( — * —  without  preconditioning 
— + —  with  the  wavelet  preconditioner). 
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1  -  Introduction 

Since  the  introduction  of  the  RWG  triangular  patch  basis  function  [1],  it  has  become  the  most 
widely  used  basis  for  solving  electromagnetic  scattering  problems  with  the  method  of  moments.  It 
has  the  highly  desirable  property  of  being  free  of  line  charges  while  being  able  to  model  arbitrarily 
shaped  surfaces.  Some  authors  have  improved  the  original  formulation  by  using  curvilinear 
triangular  patches  [2]  and  incorporating  edge  conditions  [3J.  As  pointed  out  in  [4],  it  is  well 
known  that  the  results  obtained  from  the  RWG  basis  arc  sometimes  very  dependent  on  the  meshing 
scheme  used.  Furthermore,  the  resulting  current  contains  discontinuities  that  prevent  its  use  for 
near  field  prediction.  Although  some  filtering  can  be  used  before  its  use  [5],  the  extra  filtering 
operation  results  in  loss  of  resolution.  These  problems  arc  due  to  the  fact  that  the  RWG  basis 
function  cannot  represent  an  arbitrary  linear  current  distribution  on  each  triangle. 

In  this  paper  we  present  a  new  set  of  triangular  patch  basis  functions.  It  retains  the  same 
properties  as  the  traditional  RWG  basis,  but  is  capable  of  representing  any  linear  current 
distribution  over  each  triangle.  With  this  improvement,  the  proposed  basis  is  much  less  sensitive 
to  the  meshing  scheme.  Moreover,  it  provides  a  better  representation  of  the  actual  current 
distribution,  leading  to  a  faster  convergence  of  the  moment  method  solution.  All  of  this  can  be 
achieved  with  about  the  same  computational  complexity  of  the  RWG  basis. 

2  -  Description  of  the  Problem 

If  a  closed  surface  is  represented  using  NF  discretized  triangular  patches,  a  set  of 
NE=3*NF/2  edges  connecting  those  triangles  is  obtained.  If  we  want  to  generate  a  set  of  basis 
functions  that  can  represent  an  arbitrary  linear  current  distribution  over  each  of  those  triangles,  we 
have,  for  each  triangle  /,  a  current  distribution  given  by: 

4(w,)=(a,,*,  +c„K  +(a„*,  +*„>'/  +c*)“y<  (i) 

where  x,  and  y,  arc  local  coordinates  in  the  directions  tangential  to  the  triangle  /.  Since  we  have  NF 
triangles,  the  number  of  degrees  of  freedom  for  this  problem  is  equal  to  6  NF. 

For  electromagnetic  problems  we  are  interested  only  in  current  distributions  whose 
components  normal  to  the  edges  are  continuous  to  avoid  the  presence  of  artificial  line  charges.  By 
imposing  this  additional  condition  for  each  edge  of  the  surface,  we  obtain  two  constraints  for  each 
edge,  thereby  reducing  the  degrees  of  freedom  from  6  NF  to  3  NF.  We  need,  consequently,  a  set 
of  3  NF  divergence-conforming  basis  functions  [6]  to  properly  solve  this  problem. 

In  [1  ],  Rao  et  al.  proposed  a  set  of  basis  functions  (RWG  basis  functions)  having  one  function 
associated  with  each  edge,  leading  to  a  set  of  only  3*NF/2  basis  functions.  Clearly,  such  set  cannot 
represent  an  arbitrary  linear  current  distribution  as  described  earlier,  meaning  that  it  is  not 
complete  in  that  sense.  This  has  already  been  pointed  out  in  [6],  but  the  solution  proposed  was  to 
increase  the  order  of  the  basis  using  quadratic  functions. 
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To  solve  this  problem  we  have  developed  a  new  set  of  linear  basis  functions  (which  we  shall 
term  the  TL  basis  functions)  where  we  have  two  vector  basis  functions  associated  with  each  edge 
A-c,  as  seen  in  Fig.  I  : 

p*  r  in  T* 

fi(r) 

p*  r  in  T~ 

where  p  55  r  —  F , ,  p  =  F  _  —  F  and  pA  c  =  rht.  —  ra .  Each  of  these  basis  functions  is 

parallel  to  one  edge  and  has  a  magnitude  linearly  proportional  to  the  distance  to  the  other  edge,  as 
shown  in  Fig.  2. 

The  TL  basis  functions  share  the  following  properties  with  the  RWG  functions: 

•  The  current  has  no  component  normal  to  the  boundary  of  the  surface  formed  by  triangles  T 
and  T; 

•  The  component  of  current  normal  to  the  edge  between  the  two  triangles  is  continuous  across 
the  edge; 

•  The  surface  divergence  of  these  functions,  which  is  proportional  to  the  surface  charge  density, 
is  constant  over  each  triangle,  having  the  same  magnitude  on  both,  but  opposite  signs. 

In  fact,  it  can  be  shown  that  the  RWG  basis  functions  can  be  obtained  by  adding  the  two  TL 
basis  functions  for  the  same  edge.  Because  of  the  properties  described,  it  will  be  shown  that  the 
same  electric  field  integral  equation  (EFIE)  formulation  used  with  the  RWG  functions  [I]  can  still 
be  used,  with  only  minor  modifications  for  the  new  TL  functions. 

As  an  example  of  the  improvement  that  can  be  obtained  with  the  TL  functions,  Fig.  3  shows 
the  best  representation  (in  the  mean  square  sense)  of  a  sinusoidal  current  on  a  rectangular  plate. 
We  used  the  RWG  and  the  TL  basis  functions  with  the  same  number  of  triangular  patches.  It  can 
be  seen  that  the  representation  obtained  from  the  TL  basis  is  significantly  better  than  that  from  the 
RWG  basis. 


=  >TT|rx«lp;  finr  (2) 

|p‘  *P.~|p;  r'mT~ 


3  ~  EFIE  Formulation 

Using  the  standard  EFIE  formulation  with  the  Galerkin  method  [I],  we  obtain 

(^'.7^,)  =  7^.7,,)+^./,,).  (3) 

where  =  ja  bdS  and  fm^  ^  denotes  the  basis  functions  in  (2)  at  the  m-th  edge.  As  in  [1], 

s 

because  of  the  properties  of  fm  at  the  edges,  the  last  term  in  (3)  can  be  rewritten  as 

(vo.4  _ )  =  V  ■  /; ,  j*ds + v  .f-t  \<s>ds  (4) 

One  difference  should  be  noted  here.  In  [1]  the  values  of  the  terms  involving  the  vector  potential 
and  the  incident  field  in  (3)  were  approximated  by  their  values  at  the  centroid  of  the  triangles.  With 
this  new  set  of  basis/testing  functions,  one  should  use  a  different  and  more  precise  approximation: 

(5) 

T  3  ,3,  'Mr.r, 


; 
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with 

rt  =  2/3 r.  + 1/6 r„  + 1/6 rc;  r2  =  l/6r#  +  2/3 r*  +  l/6rc;  r,  =  l/6r,  + 1 /6r,  +  2/3rc. 

This  is  necessary  because  we  now  have  3  degrees  of  freedom  for  each  triangle,  so  we  need  at  least  l 

3  testing  points.  This  approximation  leads  to  exact  results  when  the  vector  potential  or  the  incident 
field  vanes  linearly. 

The  remaining  steps  to  derive  the  matrix  equation  are  similar  to  those  in  [1J.  The  integrals 
related  to  the  scalar  and  vector  potential  determination,  after  transforming  to  a  local  system  of 
normahred  area  coordinates,  are  the  same  as  the  RWG  formulation,  so  the  modifications  needed  to 
convert  from  the  RWG  basis  to  the  TL  basis  are  fairly  straightforward. 

4  -  Numerical  Results 

To  compare  the  new  TL  basis  against  the  RWG  basis,  we  have  analyzed  the  scattering  by  a 
square  conducting  plate.  Fig.  4  and  Fig.  5  show  respectively  the  predicted  RCS  and  the  cross 
polarization  level  for  a  square  plate  of  side  IX  at  normal  incidence,  plotted  as  a  function  of  the 
number  of  basis  functions  used  (3*NF/2  for  RWG  and  3  NF  for  TL).  We  can  see  a  much  faster 
convergence  rate  with  the  TL  basis  than  with  the  RWG  basis.  Fig.  6  shows  the  current 
distributions  obtained  using  a  discretization  of  72  triangular  patches  for  the  two  types  of  basis 
functions.  We  see  veiy  abrupt  variations  of  the  current  along  the  transverse  direction  for  the  RWG 
basis.  For  the  TL  basis,  on  the  other  hand,  the  current  behavior  is  smoother,  as  expected.  It  also 
reaches  higher  values  close  to  the  side  edges. 

Other  examples,  such  as  scattering  by  a  sphere,  have  also  been  analyzed  and  will  be  shown 
during  the  presentation.  In  all  the  examples  we  have  analyzed,  we  always  obtain,  for  the  same 
meshing  scheme,  a  better  result  using  the  TL  basis  versus  the  RWG  basis  in  terms  of  RCS 
prediction.  In  terms  of  current  distribution,  the  comparative  performance  is  even  better  for  the  TL 
basis.  Therefore  this  new  set  of  basis  functions  represents  an  improvement  over  the  RWG  basis. 
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Fig.  I  Triangle  pair  and  geometrical  parameters  Fig.  2  New  improved  triangular  patch  basis 
aSSOC,ated*  function  (of  type  I )  intensity  and  direction. 


Fig.  3  (a)  Sinusoidal  current  representation  using  the  new  TL  basis  functions  and  (b)  the  original 
WG  basis  functions  over  a  square  plate,  using  72  triangular  patches.  Figure  (c)  shows  the  A-A’ 
cut  of  both  representations. 


no.  of  basis  functions 

Fig.  4  Predicted  RCS  for  normal  incidence  at  a 
square  plate  of  side  1 X,  as  a  function  of  the 
number  of  basis  functions  used. 
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Fig.  5  Predicted  cross-polarization  for  normal 
incidence  at  a  square  plate  of  side  iX,  as  a 
function  of  the  number  of  basis  functions  used. 


(a)  (b) 

Fig.  6  Current  distribution  over  a  PEC  square  plate  obtained  using  a  72  triangular  patch 
discretization  for  RWG  basis  functions  (a)  and  TL  basis  functions  (b). 
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I.  Introduction 

It  is  well  known  that  the  mutual  coupling  between  antenna  elements  in  an 
array  can  strongly  affect  Jheir  radiation/receiving  characteristics  [1,2].  For 
direction  finding  applications,  the  result  can  be  very  sensitive  to  mutual  coupling 
and  such  effect  needs  to  be  properly  accounted  for  [3],  A  common  way  to 
describe  the  coupling  effect  in  the  array  signal  processing  community  is  through 
the  use  of  a  coupling  matrix.  It  relates  the  active  element  patterns  of  the  individual 
elements  m  the  presence  of  the  array  environment  to  the  idealized,  free-standing 
element  patterns.  For  simple  antenna  structures  such  as  an  array  of  dipoles  the 
coupling  matrix  is  a  very  efficient  description  of  the  coupling  effect.  However,  for 
arrays  with  more  complex  elements  or  for  arrays  mounted  on  complex  platforms 
such  a  simple  model  is  no  longer  adequate. 

In  this  paper,  we  propose  an  improved  model  to  describe  the  active  element 
patterns  of  an  array  that  includes  both  the  mutual  coupling  and  platform  effects. 
The  key  idea  is  to  represent  the  additional  scattering  effect  from  the  platform  by 
point  scatterers  [4J.  Based  on  the  proposed  model,  we  next  devise  an  algorithm  to 
extract  the  model  parameters  from  the  observed  array  response  (i.e.,  the  active 
element  patterns)  at  a  limited  number  of  angles.  The  algorithm  entails  using  a 
matching  pursuit  approach  [5]  to  iteratively  extract  the  model  parameters.  At  each 
stage  of  the  iterauon,  the  amplitude  and  position  of  the  next  strongest  scatterer  is 
determined.  The  process  continues  until  the  difference  between  the  model  and  the 
actual  response  becomes  negligible.  Computer  simulation  results  show  that  the 
actual  array  response  can  be  accurately  recovered  from  the  improved  model. 
Consequently,  this  model  can  be  used  to  effectively  interpolate  complex  array 
responses  from  either  measurement  or  computation  made  at  a  limited  number  of 


II.  Formulation 

The  mutual  coupling  effect  in  an  antenna  array  is  commonly  modeled  by  [3]: 

Afrw  CA^  (1) 

where  A,™  is  the  actual  array  response  matrix,  A^  is  the  ideal  array  response 
matrix  in  the  absence  of  mutual  coupling,  and  C  is  an  angular-independent 
coupling  matrix.  The  rows  of  Am*  are  the  active  element  patterns  of  the  array 
elements  versus  angle.  Since  A^  is  known  based  on  the  array  geometry,  if 
is  given,  then  C  can  be  determined  by  using  the  pseudo-inverse  concept  as: 
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C  =  A,„AL(A,^aL)-'  (2) 

For  simple  antennas  whose  current  distributions  have  the  same  shape  but  different 
amplitudes,  C  can  be  found  from  at  only  a  few  observation  angles.  Once  C  is 
known,  can  be  interpolated  to  very  fine  angular  granularity  based  on  the 
model. 

When  a  platform  is  present,  the  above  simple  model  is  no  longer  valid.  We 
have  previously  proposed  an  effective  point  scatterer  model  to  describe  the 
scattering  off  a  complex  platform  [4J.  Using  this  concept,  the  scattering  from  a 
platform  can  be  equivalently  replaced  by  the  scattering  from  a  set  of  point 
scatterers.  Thus  our  proposed  model  for  the  actual  array  response  including  both 
the  mutual  coupling  and  platform  effects  can  be  written  as 

A w  =  C(A^  +  PA,to )  =  [C  Cpj^’j  (3) 

where  Apu,  is  the  total  phase  delay  due  to  the  different  point  scatterers  and  P  is  an 
amplitude  matrix  representing  the  scattering  strength  from  each  point  scatterer  to 
each  antenna  element.  In  this  model,  the  angular  dependency  due  to  the  platform 
effect  is  accounted  for  by  the  position  of  the  point  scatterers.  If  we  know  the 
position  of  the  point  scatterers,  the  matrix  C  and  P  can  be  easily  determined  by 
the  same  pseudo-inverse  procedure  as  in  (2). 

To  locate  the  position  of  the  point  scatterers,  we  employ  the  matching  pursuit 
algorithm  [5],  During  each  stage  of  the  process,  we  search  in  2D  xy-space  for  the 
best  position  of  the  next  point  scatterer.  This  is  accomplished  by  assuming  a 
point  scatterer  position,  incorporating  the  associated  phase  delay  into  the  last  row 
of  Aput  in  (3)  and  solving  for  the  corresponding  C  and  P  matrices.  The  location 
which  results  in  the  minimum 

lA^-A^u.y)!3  (4) 

is  then  chosen  as  the  next  point  scatterer.  This  process  is  similar  to  the  back 
projection  procedure  [6],  except  that  only  the  best  point  scatterer  is  retained.  This 
process  is  then  repeated  as  many  times  as  necessary  until  the  remainder  energy  is 
small  enough  to  be  neglected.  At  the  end  of  the  iteration  process,  both  the 
positions  of  the  point  scatterers  and  the  coefficients  matrices  C  and  P  are 
determined. 

III.  Results 

To  validate  the  model,  we  simulate  the  array  response  for  the  structure  shown 
in  Fig.  1  using  NEC.  The  array  consists  of  4  half-wave  dipoles  with  an  element 
separation  of  0.4  wavelength.  The  platform  is  a  plate  approximated  by  20  closely 
placed  thin  wires.  The  actual  array  response  is  computed  at  the  incident  angles 
from  0=  -90°  to  90°.  We  take  the  response  from  -90°  to  90°  at  a  step  of  10°  to 
perform  the  modeling  process  described  in  the  previous  section.  The  first  5 
extracted  point  scatterers  are  plotted  in  Fig.  2.  Over  99%  of  the  energy  of  the 
array  response  is  modeled  with  only  5  point  scatterers  in  this  case.  Fig.  3(a)  shows 
the  actual  array  response  computed  by  NEC  at  an  angular  granularity  of  1°.  Fig. 
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3(b)  is  the  array  response  predicted  by  the  model.  The  agreement  between  the  two 
is  very  good.  As  a  comparison,  we  compute  the  coupling  matrix  directly  from  (2) 
using  the  computed  result  at  all  incident  angles  (angular  granularity  of  1°).  The 
array  response  resulting  from  this  model  is  plotted  in  Fig.  3(c).  Even  though  the 
platform  is  quite  small  and  many  more  incident  angles  are  used,  the  predicted 
response  misses  some  major  characteristics.  In  direction  finding  applications,  the 
difference  can  result  in  significant  errors  in  the  determining  the  direction  of 
arrival. 

In  the  next  example,  we  look  at  a  more  complex  platform  as  shown  in  Fig.  4. 
Since  the  platform  is  10  wavelengths  in  length,  we  employ  FISC  [7],  which  is  a 
method  of  moment  (MoM)  code  based  on  the  fast  multipole  method  (FMM),  to 
simulate  the  problem.  We  compute  the  actual  array  response  from  -180°  to  180°. 
The  point  scatterers  are  extracted  using  the  response  at  the  interval  of  10°.  The 
first  15  point  scatterers  are  plotted  in  Fig.  5.  More  than  98%  of  the  energy  is 
covered  by  the  model.  The  actual  array  response  and  the  model  predicted 
response  at  an  angular  granularity  of  1°  are  plotted  in  Fig.  6(a)  and  6(b), 
respectively.  Again  the  agreement  is  very  good,  demonstrating  that  the  proposed 
model  is  effective  in  modeling  the  active  element  patterns  on  a  complex  platform. 

IV.  Conclusions 

In  this  paper,  a  model  is  proposed  to  characterize  both  the  array  mutual 
coupling  and  platform  effects.  The  platform  is  modeled  by  point  scatterers  and 
their  positions  are  determined  using  a  matching  pursuit  algorithm.  The 
effectiveness  of  the  model  is  demonstrated  by  numerical  results.  Using  this 
model,  we  can  describe  the  complex  array  response  from  measurement  or 
simulation  data  collected  at  a  limited  number  of  incident  angles. 
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(b)  top  view 

Fig.  4.  Array  -  platform  geometry  in 
example  2 


Fig.  5.  Extracted  point  scatterers 
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L  Introduction 

It  has  recently  been  demonstrated  that  the  pre-defined  wavelet  packet  (PWP) 
basis  is  very  efficient  for  representing  electromagnetic  integral  equation  [1].  The 
PWP  basis  is  designed  from  a  wavelet  packet  decomposition  tree  that  grows 
along  the  free-space  wave  number  ko.  It  is  found  that  the  non-zero  elements  in 
the  PWP-transformed  matrix  grow  at  a  rate  of  about  0(NlJ)  for  small  problem 
sizes,  and  the  rate  tends  to  0(NlogN)  for  large  problem  sizes.  Since  the  PWP- 
represented  moment  matrices  are  very  sparse  and  diagonal-dominated,  an 
efficient  preconditioner  can  be  constructed  from  the  PWP-  transformed  moment 
matrix  for  iterative  solution  of  moment  equations.  The  preconditioning  is 
especially  important  for  near-resonant  scatterers  such  as  cavity  structures,  as  the 
convergence  rate  of  solving  these  problems  using  iterative  solvers  could  be  very 
slow  [2]. 

In  this  work  we  aim  to  construct  an  effective  preconditioner  for  moment 
equations  from  the  PWP-transformed  moment  matrix.  The  preconditioning  can 
be  implemented  in  either  the  transform  domain  or  the  regular  space  domain. 
However,  the  computational  complexity  and  memory  requirement  are  at  least  on 
the  order  of  0(N2)  to  transform  the  moment  equation  from  the  space  domain  to 
the  wavelet  packet  domain  [1],  Therefore,  we  devise  a  scheme  to  compute  the 
PWP  preconditioner  directly  from  the  PWP  basis  functions  and  implement  the 
preconditioning  operation  in  the  space  domain.  Numerical  results  show  that  the 
PWP  preconditioner  is  effective  in  accelerating  the  convergence  rate  of  iterative 
solution  for  deep  cavity  structures.  We  also  demonstrate  that  the  total 
computational  complexity  and  memory  costs  for  the  preconditioner  can  be  kept 
to  0(NIogN),  making  it  compatible  with  fast  matrix-vector  multiplication 
algorithms  such  as  the  multi-level  fast  multiple  method  (MLFMM)  [3]. 

2.  PWP-Based  Approximate  Inverse  Preconditioner 

We  consider  a  left-preconditioned  system  as  follows: 

[P][Z]J  =  [P]E  (1) 

where  [Z]  is  the  space-domain  moment  matrix,  J  is  the  unknown  induced 
current  vector,  E  is  the  excitation  vector,  and  [P]  is  the  preconditioner  for  the 
moment  equation.  To  achieve  the  best  preconditioning  effect,  [P]  should  be  as 
close  to  the  inverse  of  [Z]  as  possible.  In  this  work  we  try  to  find  such  a 
preconditioner  based  on  the  PWP-transformed  moment  matrix  [Z].  It  can  be 
represented  as: 
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[Z]  =  [Mf[Z][MJ  (2) 

where  [M]T  denotes  the  transpose  of  [M]  and  is  the  PWP  transformation  matrix. 
From  (2),  noting  that  [M]  is  a  unitary  orthogonal  matrix,  we  have: 

[Z]'  =  (M][Z]'[M]t  (3) 

Since  the  PWP-transformed  moment  matrix  [Z]  is  very  sparse  and  diagonal 
dominant,  we  approximate  [  Z  ]  by  a  block-diagonal  matrix  [  Z  bd]  that  consists 
of  the  near-diagonal  terms  of  [  Z  ].  The  approximate-inverse  preconditioner  [P] 

defined  in  (1)  is  then  constructed  by  inverting  the  block-diagonal  matrix  [ZM] 
and  transforming  the  resulting  matrix  back  to  the  space  domain: 

[P]  =  [M][Zbdrl[M]T  (4) 

Therefore,  the  preconditioned  moment  equation  becomes: 

[Ml  IZbd]'1  [M]T[Z]  J  =  [M]  [ZM l'1  [M]tE  (5) 

The  inverse  of  the  block  diagonal  matrix  [ZM]  is  simply  the  inverse  of  its 
individual  diagonal  blocks,  and  the  inverted  matrix  [Zm]'1  remains  block 
diagonal.  By  properly  choosing  the  block  sizes  of  [ZM],  we  can  limit  the 
computational  cost  of  the  inversion  while  maintaining  the  effectiveness  of  the 
preconditioner.  Our  approach  is  to  divide  the  diagonal  area  of  [Z]  into  three 
parts  (Fig.  1):  (i)  a  block  centered  about  the  spectral  frequency  ko  with  block  size 
Mt,  (ii)  a  set  of  diagonal  blocks  in  the  remaining  upper-left  region  with  block 
size  M2,  and  (iii)  a  set  of  diagonal  blocks  in  the  lower-right  region  with  block 
size  M3.  Because  most  of  the  energy  in  [Z]  is  concentrated  around  the  ko 
spectral  frequency,  we  let  Mi  grow  a  rate  of  */n  with  the  problem  size  N.  Next, 
with  the  chosen  pattern  of  [  Z  bdl,  we  compute  the  non-zero  elements  one-by-one 
using  the  PWP  basis  functions  directly  [1].  Since  the  number  of  non-zero 
elements  in  [  Zw]  is  about  0(N)  and  the  average  cost  to  compute  an  element  is 
about  O(logN),  we  can  show  that  the  total  computational  cost  and  memory 
requirement  to  generate  and  invert  [ZM]  are  under  0(NlogN)  [4],  Finally,  the 
PWP  transformation  matrix  [M]  in  (5)  is  sparse  with  about  0(NlogN)  non-zero 
elements.  Consequently,  the  computation  complexity  and  memory  requirement 
for  the  PWP  preconditioning  operation  can  be  kept  to  0(NlogN). 

3.  Numerical  Results 

To  test  the  performance  of  the  PWP  preconditioner,  we  use  an  open-ended 
inlet  shown  in  Fig.  2(a)  as  the  scatterer.  The  combined  field  integral  equation 
(CFIE)  is  used  to  construct  the  moment  equation.  Fig.  2(b)  is  the  moment 
matrix  [  Z]  represented  by  the  PWP  basis  for  the  inlet  structure  with  dimensions 
/:m:f=15:l:l  and  N=256.  We  observe  that  the  matrix  is  sparse  and  diagonal- 
dominant,  with  the  strongest  elements  located  around  ko  and  the  low  spectral 
frequency  areas. 

Fig.  3  shows  the  iteration  number  required  to  solve  the  PWP  preconditioned 
system  versus  the  problem  size  with  /:m:f=15:l:l.  In  the  PWP  preconditioner. 
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we  choose  M2=32,  M3=l,  and  let  M|  grow  at  a  rate  of  V#  starting  from  32  at 
N=128.  The  problem  size  N  is  increased  by  proportionally  increasing  the 
scatterer  size.  The  iterative  solver  used  is  the  conjugate-gradient  squared  (CGS) 
method,  and  the  stopping  criterion  is  when  the  relative  residual  error  is  less  than 
10"*.  For  comparison,  we  also  show  the  results  from  a  simple  block-diagonal 
preconditioner  in  the  space  domain.  We  find  that  the  iteration  number  grows 
rapidly  with  the  problem  size  if  no  preconditioning  is  applied.  With  the  PWP 
preconditioner  applied,  the  iteration  numbers  are  significantly  reduced  and  they 
increase  more  slowly  with  problem  size.  Furthermore,  the  PWP  preconditioner 
performs  better  than  the  equivalent  space  domain  preconditioner. 

Fig.  4  shows  the  convergence  behavior  of  the  preconditioned  system  for  a 
fixed  problem  size  (N=512)  with  the  BiCGSTAB  algorithm  as  the  iterative 
solver.  The  residual  error  is  plotted  as  a  function  of  the  number  of  iterations  for 
the  moment  equation  without  any  preconditioning,  or  with  preconditioning  in 
the  space  domain  or  the  PWP  basis  domain.  We  observe  that,  with  the  PWP 
preconditioning,  the  relative  residue  decreases  the  fastest.  In  addition,  the 
convergence  behavior  is  more  stable.  Similar  convergence  behaviors  are  found 
for  larger  problem  sizes  and  for  other  near-resonant  configurations. 

4.  Conclusions 

Due  to  the  vanishing  moments  of  wavelet  basis  functions,  the  PWP-based 
moment  matrix  is  sparse  and  diagonally  concentrated.  Consequently,  an 
effective  preconditioner  can  be  more  easily  constructed  than  that  based  on  the 
space-domain  moment  matrix.  Numerical  results  showed  that  the  iteration 
numbers  for  the  PWP-preconditioned  moment  equations  are  significantly 
smaller  and  grow  at  a  lower  rate  than  those  without  preconditioning  or 
preconditioned  using  a  space-domain  preconditioner.  The  complexity  for  the 
construction  and  preconditioning  operation  of  the  PWP  preconditioner  is  kept 
under  O(NlogN)  in  both  computational  cost  and  memory  requirement. 
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Figure  1.  The  block  structure  used  to  construct  [  Z  w] 


by  PWP  basis[  Z  ],  for  the  inlet  with  N=256  (in  logarithmic  scale). 


Figure  3.  Iteration  numbers  vs.  problem  Fig.4.  Convergence  behaviors  of  the 
sizes  for  solving  moment  equations  various  preconditioned  systems 

using  different  preconditioning  methods  versus  iteration  numbers  for  PWP 

and  other  preconditioning  methods 
with  N=512. 
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